A gravity anomaly and vertical gradient joint inversion method, system and electronic equipment

Abstract

The application discloses a gravity anomaly and vertical gradient joint inversion method and system and electronic equipment, and belongs to the cross technical field of geophysical exploration and deep learning. The application stacks gravity anomaly and gravity vertical gradient data along the channel dimension into a double-channel input tensor, inputs the double-channel input tensor into an UNet++ network with the first layer convolution input channel number of the encoder being 2, makes the network automatically learn the complementary feature fusion weight of the two kinds of data, realizes joint inversion of deep and shallow collaborative constraints, adopts a MSE and MAE weighted combination loss function to balance global shape recovery and local detail optimization, independently normalizes the double channels respectively, and persists the parameters, so that the gradient leading problem caused by the data magnitude difference is avoided, and a double-core matrix is loaded during reasoning to realize joint evaluation. The application solves the problems of difficult weight selection, large calculation amount and insufficient single data source inversion precision in the traditional joint inversion, and the precision of the application is significantly better than that of the single data source inversion in the deep and shallow parts.

Description

Technical Field

[0001] This invention relates to the field of geophysical exploration and deep learning interdisciplinary technology, and in particular to a method, system and electronic device for joint inversion of gravity anomalies and vertical gradients. Background Technology

[0002] Gravity exploration is a crucial method in geophysical exploration. By measuring gravity anomalies at or near the Earth's surface, it infers subsurface density distribution and is widely used in oil and gas exploration, mineral surveys, engineering geological investigations, and crustal structure research. In gravity exploration, gravity anomaly Δg and vertical gravity gradient are key components. (VGG) are two different physical observations that complement each other in their response characteristics to subsurface density volumes. Gravity anomaly Δg is sensitive to large-scale density volumes at depth, reflecting the overall mass effect of the density volume, but its spatial resolution is limited, making it difficult to distinguish between superimposed anomalies at depth and shallow depth. The gravity kernel matrix G contains ln and arctan terms; the ln term provides a far-field contribution, enhancing the response of deep elements, and the decay of kernel matrix elements with depth z is approximately O(1 / z²). Gravity vertical gradient (VGG) is sensitive to small-scale density volumes at shallow depth, has high vertical resolution, and can better distinguish vertically superimposed density volumes, but its response to depth is weak. The VGG kernel matrix G_vgg contains only an arctan term, decays faster in the far field, and the decay of kernel matrix elements with depth z is approximately O(1 / z³). The contribution of deep elements to VGG is smaller than that of gravity anomalies, while the contribution of shallow elements is relatively larger. Their forward modeling formulas are as follows: The unit is mGal; The unit is Eötvös (E). Where G and G_vgg are the gravity kernel matrix and VGG kernel matrix, respectively, and ρ is the density contrast (kg / m³). 3 ), where m is the binary model vector. The two types of observation data share the same density model (same source), differing only in their kernel matrices. The gravity kernel matrix G is based on the Nagy rectangular gravity anomaly analytical formula (containing ln and arctan terms), while the VGG kernel matrix G_vgg contains only the arctan term (partial derivative of gravity anomaly in the z-direction).

[0003] The core idea of ​​joint inversion is to simultaneously use the two types of observation data to constrain the same density model, with the deep constraint derived from Δg and the shallow constraint from VGG, achieving high-precision inversion across the entire depth range. Existing joint inversion techniques mainly include the following categories: 1. Traditional weighted joint inversion method: This method uses a weighted objective function to simultaneously fit two types of observation data. For example, the joint inversion method proposed by Gallardo-Delgado et al. (2003) is a typical representative of this field. The specific function is as follows:

[0004] Where W_d is the data weight matrix, and R(m) is the regularization term. and These are data weight parameters. This method has the following drawbacks: (1) Difficulty in selecting data weights: , The selection of the data source significantly impacts the inversion results, with optimal weights varying greatly across different scenarios, and an adaptive selection method is lacking. If the weight is too large, the data source dominates the inversion results; if the weight is too small, the contribution of the data source is ignored, requiring repeated experimentation to determine the optimal weight. Typically, , The search requires 10 to 20 inversion experiments, which is computationally expensive; (2) High computational cost: Joint inversion requires processing two sets of observation data simultaneously, doubling the scale of iterative solutions and resulting in low computational efficiency. Joint inversion of a 32×32×16 grid takes approximately 10-30 minutes per iteration, which is 2-3 times that of single-source inversion; (3) Regularization parameter coupling: data weights , Coupled with the regularization parameter λ, the parameter space search is difficult and requires a three-dimensional parameter space. Search for the optimal combination.

[0005] 2. Sequential Inversion Method: This method first inverts one type of data, then uses the result as the initial model or constraint for inverting another type of data. For example, the cross-gradient joint inversion proposed by Fregoso & Gallardo (2009) belongs to this category. This method has the following drawbacks: (1) One-way information transmission: Errors in the initial inversion results will be transmitted to subsequent inversions, making true two-way constraints impossible. For example, errors in the deep model obtained by inverting Δg first will affect the shallow results of VGG inversion; (2) Dependence on inversion order: Different orders lead to different results and lack objectivity. Inverting Δg first focuses on deep accuracy, while inverting VGG first focuses on shallow accuracy. The results of the two orders may differ significantly.

[0006] 3. Deep learning-based single-data-source inversion methods: Existing deep learning inversion methods only use a single data source (gravity anomaly or VGG), failing to fully utilize the collaborative constraints of multiple data sources. These methods have the following drawbacks: (1) Insufficient depth resolution: Inversion using only a single data source (gravity anomaly or VGG) does not fully utilize the collaborative constraints of multiple data sources. When inversion using only VGG, the recovery of the deep density volume is poor (VGG has a weak response to the deep part, and the deep R² decreases by about 15%). (2) Shallow boundary blurring: When only Δg is used for inversion, the shallow boundary is blurred (Δg has low spatial resolution, and the shallow R² decreases by about 10%). (3) Lack of multi-data fusion mechanism: The dual-channel input architecture was not designed, making it impossible to process two types of physical field data simultaneously. Even if two networks are trained separately and then the results are fused, end-to-end multi-data collaborative optimization cannot be achieved. Summary of the Invention

[0007] The purpose of this invention is to overcome the technical problems existing in the prior art and to provide a method, system and electronic device for joint inversion of gravity anomalies and vertical gradients.

[0008] The objective of this invention is achieved through the following technical solution: A first aspect of the present invention provides a method for joint inversion of gravity anomalies and vertical gradients, comprising the following steps: S1. Discretize the underground space into multiple cubic elements and calculate the gravity kernel matrix and gravity vertical gradient kernel matrix respectively; S2. Generate a density model, and perform forward modeling using the gravity kernel matrix and the gravity vertical gradient kernel matrix respectively to obtain gravity anomaly data and gravity vertical gradient data corresponding to the same density model; S3. Construct a dual-channel input UNet++ network to stack gravity anomaly data and gravity vertical gradient data along the channel dimension into a dual-channel input tensor; S4. The dual-channel input UNet++ network is trained using a joint inversion loss function; S5. Input the gravity anomaly data and gravity vertical gradient data to be inverted into the trained dual-channel input UNet++ dense skip connection network to obtain the prediction density model.

[0009] In some embodiments, the dual-channel input UNet++ network adopts a dense skip connection structure, in which the encoder downsamples layer by layer to extract multi-scale features, and the decoder fuses the features of each layer through dense skip connections; the number of input channels of the first convolutional layer of the encoder is 2, which is used to simultaneously receive gravity anomaly data and gravity vertical gradient data.

[0010] In some embodiments, the gravity kernel matrix is ​​calculated based on the Nagy rectangular gravity anomaly analytical formula, and the kernel function includes ln terms and arctan terms; the gravity vertical gradient kernel matrix is ​​obtained by taking the z-direction partial derivative of the gravity formula, and the kernel function contains only arctan terms.

[0011] In some embodiments, the shape of the dual-channel input tensor is (2, H, W), with one channel inputting gravity anomaly data and the other channel inputting gravity vertical gradient data; the size of the first layer convolutional kernel of the encoder is (2, 3, 3), and the kernel weights are automatically learned through backpropagation to obtain the optimal fusion weights for gravity anomaly data and gravity vertical gradient data.

[0012] In some embodiments, the joint inversion loss function is a weighted combination of mean squared error (MSE) and mean absolute error (MAE): L = α·MSE + β·MAE, where α and β are weighting coefficients, the mean square error (MSE) is used to eliminate large deviations to restore the global shape, and the mean absolute error (MAE) is used to reduce small deviations to optimize local details.

[0013] In some embodiments, during training in step S4, the z-score of the dual-channel input is independently normalized and the normalization parameters of each channel are persistently saved.

[0014] In some embodiments, gravity anomaly data and gravity vertical gradient data corresponding to the same density model are paired and assembled according to their numbers to generate joint data, ensuring that the two types of observation data share the same field source model.

[0015] In some embodiments, after obtaining the predicted density model, the gravity kernel matrix and the gravity vertical gradient kernel matrix are loaded, and the predicted gravity anomaly and the predicted gravity vertical gradient are calculated based on the forward model of the predicted density model. They are then compared and verified with the observation data to achieve joint evaluation of the dual kernel matrix.

[0016] A second aspect of the present invention provides a joint inversion system for gravity anomalies and vertical gradients, comprising: The kernel matrix generation module is used to discretize the underground space into multiple cubic units and calculate the gravity kernel matrix and gravity vertical gradient kernel matrix respectively. The forward modeling module is used to generate density models. It performs forward modeling calculations using the gravity kernel matrix and the gravity vertical gradient kernel matrix to obtain gravity anomaly data and gravity vertical gradient data corresponding to the same density model. The dual-channel input UNet++ network building module is used to build a dual-channel input UNet++ network, which stacks gravity anomaly data and gravity vertical gradient data along the channel dimension into a dual-channel input tensor; A dual-channel input UNet++ network training module is used to train the dual-channel input UNet++ network using a joint inversion loss function; The inversion module is used to input the gravity anomaly data and gravity vertical gradient data to be inverted into the trained dual-channel input UNet++ dense skip connection network to obtain the prediction density model.

[0017] A third aspect of the present invention provides an electronic device, including a memory and a processor, wherein the memory stores computer instructions executable on the processor, and the processor executes a gravity anomaly and vertical gradient joint inversion method as described in the first aspect when executing the computer instructions.

[0018] It should be further noted that the technical features corresponding to the above-mentioned options and embodiments can be combined or substituted with each other to form new technical solutions without conflict.

[0019] Compared with the prior art, the beneficial effects of the present invention are: 1. Deep-shallow co-constraints improve full-depth inversion accuracy: Dual-channel input enables the network to simultaneously acquire gravity anomalies (deep information) and VGG (shallow information). Compared to single-source inversion, the deep R² increases from 0.88 to 0.92, and the shallow R² increases from 0.82 (Δg only) or 0.90 (VGG only) to 0.94, with an overall R² of 0.93. This advantage stems from the deep-shallow co-constraint architecture of the dual-channel input—the Δg channel provides deep constraints (ln term far-field contribution), and the VGG channel provides shallow constraints (arctan term near-field sensitivity). 2. Automatic learning of data fusion weights, eliminating the need for manual parameter tuning: Traditional joint inversion requires manual setting of data weights. , (Typically requiring 10-20 trials), this method automatically learns the optimal fusion method through the first-layer convolutional kernel (2×3×3), avoiding the difficulty of weight selection. The convolutional kernel weights are automatically optimized through backpropagation, with deeper regions relying more on Δg channel features and shallower regions relying more on VGG channel features, achieving adaptive hierarchical fusion. This advantage stems from the dual-channel input and the automatic feature fusion mechanism of convolution; 3. The joint loss function balances global and local accuracy: The MSE+MAE combined loss takes into account both global morphology and local details. Compared with a single loss function, the joint inversion result has a more accurate overall morphology and sharper boundaries. MSE's quadratic penalty quickly eliminates large deviations (such as overall density volume position shifts), while MAE's linear penalty finely adjusts small deviations (such as pixel-level boundary shifts). This advantage stems from the weighted combined loss design of MSE+MAE. 4. Significantly Improved Inversion Efficiency: After training, inference requires only one forward propagation (approximately 5ms per sample, GPU environment), simultaneously outputting joint inversion results, representing a 4-5 order of magnitude improvement compared to traditional joint iterative methods (10-30 minutes / cycle). This advantage stems from the end-to-end mapping of deep learning—traditional methods require iteratively solving large linear equation systems (computational complexity O(N_cell²·N_iter)), while deep learning methods encode the inversion process into network weights (forward propagation has a fixed computational cost); 5. Fully compatible with single data source framework: Joint inversion only requires modification of the number of input channels (1→2) and the loss function (Dice+BCE→MSE+MAE). The main network architecture is completely consistent with that of single data source inversion, facilitating engineering implementation and maintenance. This advantage stems from the unified network architecture design. Attached Figure Description

[0020] Figure 1 This is a flowchart illustrating a method for joint inversion of gravity anomalies and vertical gradients according to an embodiment of the present invention; Figure 2 This is a flowchart illustrating the overall technical process of joint inversion according to an embodiment of the present invention; Figure 3 This is a diagram illustrating the dual-channel input UNet++ architecture of an embodiment of the present invention; Figure 4 This is a schematic diagram of the training / validation / test loss curves shown in an embodiment of the present invention; Figure 5 This is a schematic diagram illustrating the joint inversion results in an embodiment of the present invention. Detailed Implementation

[0021] The technical solution of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. The components of the embodiments of this application described and shown in the accompanying drawings can generally be arranged and designed in various different configurations. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0022] It should be noted that the defects in the solutions in the prior art are all the results of the inventors' practice and careful research. Therefore, the discovery process of the above problems and the solutions proposed by the embodiments of this application in the following text should be the inventors' contributions to this application in the process of invention and creation, and should not be understood as technical content known to those skilled in the art.

[0023] Based on the above statements, and since this embodiment involves gravity inversion-related technologies, the following explanations are provided to clarify the objectives, technical solutions, and advantages of this application embodiment: 1. UNet++: Nested U-Net architecture, dense skip connections enable multi-scale feature fusion.

[0024] 2. Dense skip connections: Each decoder node receives features from all preceding encoder layers, unlike the single-layer skip connections in standard U-Net.

[0025] 3. Semantic gap: The information gap between shallow features (high resolution, low semantics) and deep features (low resolution, high semantics) of the encoder.

[0026] 4. Dice Loss: A loss function based on the Dice coefficient, which is robust to class imbalance.

[0027] 5. VGG: Vertical Gravity Gradient, the z-direction derivative of the gravitational anomaly.

[0028] 6. Kernel Matrix: The discretization matrix of the forward operator, which realizes a linear mapping from the model space to the data space.

[0029] 7. DACE Kriging: A Kriging interpolation method based on computer experimental design and analysis.

[0030] 8. z-score normalization: (x-μ) / σ standardization makes the data mean 0 and the standard deviation 1.

[0031] 9. ELU: Exponential Linear Unit, with a non-zero gradient in the negative region, avoiding the ReLU dead zone.

[0032] 10. AdamW: An Adam optimizer that decouples weight decay, resulting in more stable regularization performance.

[0033] In view of the technical problems pointed out in the background art, the present invention provides the following embodiments: In one exemplary embodiment, a method for joint inversion of gravity anomalies and vertical gradients is provided, such as... Figure 1 As shown, it includes the following steps: S1. Discretize the underground space into multiple cubic elements and calculate the gravity kernel matrix and gravity vertical gradient kernel matrix respectively; S2. Generate a density model, and perform forward modeling using the gravity kernel matrix and the gravity vertical gradient kernel matrix respectively to obtain gravity anomaly data and gravity vertical gradient data corresponding to the same density model; S3. Construct a dual-channel input UNet++ network to stack gravity anomaly data and gravity vertical gradient data along the channel dimension into a dual-channel input tensor; S4. The dual-channel input UNet++ network is trained using a joint inversion loss function; S5. Input the gravity anomaly data and gravity vertical gradient data to be inverted into the trained dual-channel input UNet++ dense skip connection network to obtain the prediction density model.

[0034] The following describes the implementation process of the present invention, such as... Figure 2 As shown, the specific steps include: Step 1. Discretization of underground space and calculation of dual-kernel matrix The underground space is discretized into Nsx×Nsy×Nbz cubic elements (in this embodiment, the default is 32×32×16, with a side length of 50m). Nsx represents the number of grid cells in the x-direction, Nsy represents the number of grid cells in the y-direction, and Nbz represents the number of vertical grid cells. The following calculations are performed: (1) Gravity kernel matrix G: Based on the Nagy rectangular gravity anomaly analytical formula, the gravity response of each observation point to each cubic unit is calculated. The shape of the gravity kernel matrix G is (1024, 16384), where 1024 is the number of observation points (32×32) and 16384 is the number of grid units (32×32×16), with units of mGal / (kg / m³), and it is persistently stored in G.npy format. The kernel function of this gravity kernel matrix G contains ln terms (natural logarithm terms) and arctan terms (arctangent terms). (2) VGG kernel matrix G_vgg: calculates the z-direction partial derivative (vertical derivative) of the gravity formula. The kernel function contains only the arctan term (arctangent term). The output unit is Eötvös (E). The shape of the VGG kernel matrix is ​​(1024, 16384), and the unit is E / (kg / m³). It is persistently stored in the format G_vgg.npy.

[0035] Step 2. Training Sample Density Model Generation A density model is generated by combining a random walk model (90%) with a regular geological model (10%). The same model generation logic is used as that used in single-data source inversion to ensure a fair comparison between joint inversion and single-data source inversion.

[0036] (1) Random walk model (accounting for 90% of the total): Starting from a 2×2×2 initial cube (a seed cube composed of 8 adjacent voxels), it expands along a random axis (one of the three axes x / y / z) and a step size (±2 voxels, i.e., ±100m). The number of walk steps is randomly selected within the range of [step1, step2] (default step1=60, step2=80). It supports the generation of 1 or 2 density volumes (50% probability each), providing rich morphological diversity. The model is stored in .dat format, with each line recording (ix, iy, iz, value), where value=1.0 indicates that the unit belongs to a density volume. The design basis for 60~80 random walk steps: if the number of steps is too small (<30 steps), the density volume will be too small, the forward modeling anomaly signal will be weak, which is not conducive to network learning; if the number of steps is too large (>100 steps), the density volume will be too large, occupying too much underground space, and the class imbalance problem will be aggravated. (2) Regular geological models (accounting for 10% of the total): contain 6 typical geological structures, each generated in equal quantities: Type I Rectangular Prism: A cuboid with a length of 8~24 blocks (400~1200m) in the x direction, a width of 8~24 blocks in the y direction, a height of 3~8 blocks (150~400m) in the z direction, and a random center position; Type II Inclined Dike: An inclined plate-like body along the xz plane, with the dip angle achieved by offsetting in the z direction, simulating an intrusive dike; Type III Vertically Separated Prisms: Two vertically separated prisms, simulating layered ore bodies; Type IV syncline prism: a syncline structure that is deep in the middle and shallow on both sides, simulating a syncline ore body; Type V Parallel Vertical Prism: Two horizontally separated parallel prisms, simulating parallel veins; Type VI Fault Prism: A vertically displaced prism that simulates the ore bodies on both sides of a fault; The design principles of the six regular models are: to cover the most common geological structure types—intrusive bodies (Type I / II), layered ore bodies (Type III), fold structures (Type IV), parallel veins (Type V), and fault structures (Type VI)—to ensure that the network learns the feature representation of regular geometric shapes and to compensate for the regular boundary shapes that the random walk model cannot cover.

[0037] Step 3. Forward modeling and joint data generation from dual data sources (1) Generate gravity anomaly data: d_grav = G·(ρ·m) and VGG data: d_vgg = G_vgg·(ρ·m), where ρ is a fixed density contrast (default 1000 kg / m³). (2) Joint data assembly: Package the d_grav and d_vgg corresponding to the same density model into a joint data file (.npz format), which contains three fields: d_grav: Gravity anomaly data, shape (1, 1024), unit mGal; d_vgg: Gravity vertical gradient data, shape (1, 1024), unit E; m: density model (multiplied by density value), shape (1, 16384), unit kg / m³.

[0038] (3) Joint data generation process: Traverse the generated gravity anomaly data and VGG data, pair them one by one according to their numbers, and assemble them into a joint data file. This design ensures that the joint data and the single data source data use the exact same density model, guaranteeing the fairness of the comparative experiment. The joint data storage directories are . / joint_data / train / , . / joint_data / val / , and . / joint_data / test / .

[0039] Step 4. Dual-channel input UNet++ network construction Building upon the UNet++ dense skip connection architecture, key modifications are made to address the dual-data source characteristics of joint inversion. Gravity anomaly and vertical gravity gradient (VGG) data are stacked along the channel dimension as a (2, H, W) dual-channel input. By modifying the number of input channels in the first layer of the encoder, the network automatically learns to fuse complementary features from both datasets, achieving joint inversion with deep and shallow part co-constraints. The dual-channel input is equivalent to increasing the number of observation constraints from N_obs to 2×N_obs, enhancing the constraints and making the network output more deterministic. Based on the dual-channel input, the UNet++ dense skip connection architecture is employed, where each decoder node receives features from all preceding encoder layers, bridging the semantic gap and improving the boundary recovery accuracy of the joint inversion. The specific structure is as follows: Figure 3 As shown: (1) Input channel modification: The number of input channels in the first convolutional layer of the encoder is changed from 1 to 2, i.e., in_channels=2, so that the network can receive both gravity anomaly and VGG data simultaneously. The kernel size of the first convolutional layer is (2, 3, 3), where the 2 input channels correspond to the 2 physical fields, the 3×3 spatial kernel extracts the local spatial pattern, and the convolutional kernel weights automatically learn the best linear combination of the two types of data; (2) Input data organization: stack d_grav and d_vgg along the channel dimension to form a two-channel input tensor with the shape (2, Nsx, Nsy) i.e. (2, 32, 32), where channel 0 is the gravity anomaly and channel 1 is VGG; (3) Encoder-decoder architecture: The 5-layer UNet++ dense skip connection structure remains unchanged. The encoder downsamples layer by layer to extract multi-scale features (number of channels 64→128→256→512→1024). The decoder fuses the features of each layer through dense skip connections. The final decoder node of Level 0 receives features from 5 preceding nodes (number of input channels start_fm×5=320). start_fm is the number of initial feature maps. In this embodiment, it is taken as 64 to achieve the most sufficient multi-scale feature fusion. (4) Output layer: Double Conv(start_fm, Nbz) → 1×1 Conv2d(Nbz, Nbz, kernel=1) → Batch Normalization layer BatchNorm → Sigmoid activation function, output geometric mask ∈[0,1]^(Nbz×Nsx×Nsy), mapping the output value to the interval [0,1], consistent with the output of single data source inversion.

[0040] Step 5. Design of Joint Inversion Loss Function The joint inversion uses a weighted combined loss function of MSE+MAE: joint_loss = α·MSE(pre_y, tru_y) + β·MAE(pre_y, tru_y) Where α = 0.5, β = 0.5. Design rationale: (1) MSE (mean squared error): MSE = (1 / N)Σ(pre_y - tru_y)², this loss imposes a second penalty on large errors, which is beneficial to quickly eliminate significant biases in the prediction model (such as the overall positional shift of the density volume) and ensure the correct overall shape of the density volume. (2) MAE (Mean Absolute Error): MAE = (1 / N)Σ|pre_y - tru_y|, this loss applies a uniform linear penalty to all errors, which is beneficial for fine-tuning small biases in the prediction model (such as boundary pixel-level offsets) and improving boundary accuracy. (3) Weighted combination: MSE focuses on the global form, while MAE focuses on local details. The two complement each other to achieve precise constraints for joint inversion; (4) Reasons for not using Dice+BCE: In joint inversion, dual-channel input already provides stronger constraint information, the network output is more stable, and the class imbalance problem is alleviated. There is no need for Dice's robust protection against class imbalance, and the MSE+MAE combination is more direct and effective. Specifically, dual-channel input is equivalent to increasing the number of observation constraints from N_obs to 2×N_obs. The enhanced constraints make the network output more certain (the prediction variance is reduced), and the impact of class imbalance is diluted—even if the density volume only accounts for 5% of the voxels, the gradient signal under dual-channel constraints is still sufficient to drive the network to segment correctly. This design idea of ​​"adaptively selecting the loss function according to the strength of the input constraints" is the core consideration for the selection of the joint inversion loss function.

[0041] Step 6. Training Strategy (1) Data preprocessing: The dual-channel inputs are z-score independently normalized: the normalized gravity anomaly data d_grav_norm = (d_grav - μ_grav) / σ_grav, and the normalized gravity vertical gradient data d_vgg_norm = (d_vgg - μ_vgg) / σ_vgg; Four normalization parameters (μ_grav, σ_grav, μ_vgg, σ_vgg) are stored and loaded during inference to ensure consistency; μ_grav (mean of gravity anomaly) and σ_grav (standard deviation of gravity anomaly) are the normalization parameters for gravity anomaly data, and μ_vgg (mean of gravity vertical gradient) and σ_vgg (standard deviation of gravity vertical gradient) are the normalization parameters for gravity vertical gradient data.

[0042] The necessity of independent normalization: Gravity anomaly data typically ranges from 0 to 10 mGal, while gravity vertical gradient data typically ranges from 0 to 100 E, with a significant difference in magnitude. If a uniform normalization parameter is used instead of independent normalization, the value of the gravity vertical gradient channel is much larger than that of the gravity anomaly channel. This causes the gravity vertical gradient channel to dominate gradient updates during network training, while the contribution of the gravity anomaly channel is ignored, resulting in a degeneration of the joint inversion into an approximate single gravity vertical gradient inversion. After independent normalization, the numerical ranges of both types of data are approximately [-3, 3], with consistent magnitudes, and the contributions of the two channels to the gradient are balanced.

[0043] Furthermore, the label model density is normalized to m / ρ (binary model vector divided by density contrast) to make the output value ∈ [0,1], matching the Sigmoid activation function. The dual-channel input data is reshaped to (2, 32, 32), and the labels are reshaped to (16, 32, 32).

[0044] The dual-channel input data is reshaped to (2, 32, 32), and the labels are reshaped to (16, 32, 32). (2) Data Loading: The GravityDataset class (gravity dataset class) returns a (image, mask) tuple in joint mode, where image (input image) is a two-channel normalized data with a shape of (2, 32, 32), and mask (mask label) is a normalized geometric mask with a shape of (16, 32, 32). In joint mode, the data file contains two fields, d_grav and d_vgg, which are read by DataLoader and stacked along the channel dimension. (3) Optimizer: AdamW, learning rate lr=3×10 -4 Weight decay weight_decay=1e-4 AdamW decouples weight decay from gradient updates, resulting in more stable regularization. (4) Early stopping mechanism: No improvement in validation set loss over consecutive patience rounds (default 20) (change less than threshold = 1 × 10). -4 If the model loses the least on the validation set, training is stopped, and the model weights with the lowest loss on the validation set are retained to prevent overfitting. (5) Checkpoint continuation training: Save the model weights (.pth format), optimizer state, 4 normalization parameters, and loss history to the _ckpt.pth file (checkpoint file) for each epoch (complete traversal of the training set), and support breakpoint resumption training.

[0045] Step 7. Reasoning and Evaluation (1) Load 4 normalization parameters to normalize the dual-channel input respectively. The model forward propagation obtains the prediction mask pred_mask∈[0,1]^(16×32×32), which is multiplied by the density value ρ to recover the density model: density_model = pred_mask × ρ (unit kg / m³), or density_model = pred_mask × ρ / 1000 (unit g / cm³). (2) 3D voxel visualization: Based on fixed color mark mapping (density 0~1.0 g / cm³, one color level for every 0.1 g / cm³), ensuring color comparability between different samples; (3) Comparison of contour lines for dual anomalies: Load the dual-core matrices G and G_vgg, and calculate two predicted anomalies (d_grav_pred = G @ (pred_mask×ρ).flatten(), d_vgg_pred = G_vgg @ (pred_mask×ρ).flatten()) based on the prediction density model. Plot contour lines with the observed anomalies d_grav_obs and d_vgg_obs respectively. Physical significance of the consistency verification of dual anomaly forward modeling: If the prediction model is correct, the two anomalies it forward models should be consistent with the two observed anomalies respectively. Any deviation of any anomaly reflects the insufficiency of the inversion model. (4) Quantitative assessment: relative error The coefficient of determination R² = 1 - Σ(m_pred - m_true)² / Σ(m_true - mean(m_true))². R² is evaluated separately for deep and shallow regions. The criterion for dividing deep and shallow regions is that the first 8 layers in the z-direction (z=0~7, depth 0~400m) are considered shallow, and the last 8 layers (z=8~15, depth 400~800m) are considered deep.

[0046] It should be noted that the core principle of the dual-channel fusion in this invention is the automatic learning of data weights. The dual-channel input allows the network to simultaneously access two types of physical field data in the first layer of the encoder, and the convolutional kernel automatically learns the complementary feature fusion method of the two types of data. (1) Channel 0 (gravity anomaly) provides deep large-scale information—gravity anomalies have a strong response to deep density volumes, with gentle spatial changes, reflecting the overall quality effect of the density volume. After the first layer of convolution in the encoder, the feature map corresponding to channel 0 has gentle spatial changes (mainly low-frequency components), containing the overall position and scale information of the deep density volume; (2) Channel 1 (VGG) provides high-resolution information on shallow density volumes—VGG has a strong response to shallow density volumes, with dramatic spatial changes, reflecting the vertical structural details of the density volumes. After the first layer of convolution in the encoder, the feature map corresponding to Channel 1 undergoes dramatic spatial changes (mainly high-frequency components), containing the boundary positions and morphological details of shallow density volumes; (3) The first convolutional kernel (2×3×3) automatically learns the optimal fusion weights for the two types of data, eliminating the need for manual setting of data weights. , Specifically, the weights of the two input channels of the convolutional kernel correspond to the contribution ratios of the two physical field data, and are automatically optimized through backpropagation. During training, the gradient signal guides the convolutional kernel weights to adjust to the optimal ratio—deep regions rely more on the features of channel 0 (Δg), while shallow regions rely more on the features of channel 1 (VGG). The network automatically achieves deep and shallow layer fusion. Subsequent convolutional layers further extract complementary features from the two types of data. The deep features fuse information from both deep and shallow regions, fundamentally solving the problem of difficult data weight selection in traditional joint inversion.

[0047] Furthermore, a quantitative comparison was made between the joint inversion of the present invention and the single data source inversion. Under the experimental conditions of 32×32×16 grid, 20,000 training samples, 300 epochs, batch_size=32, and NVIDIA GPU, the accuracy improvement of the joint inversion compared with the single data source inversion is shown in Table 1.

[0048] Table 1. Quantitative Comparison Results of Joint Inversion and Single Data Source Inversion

[0049] The results show that the joint inversion of this invention improves the R² value from 0.88 (gravity anomaly only) to 0.92 (a 4.5% improvement) in the deep region and from 0.82 (gravity anomaly only) or 0.90 (gravity vertical gradient only) to 0.94 in the shallow region (a 14.6% improvement compared to gravity anomaly only, and a 4.4% improvement compared to gravity vertical gradient only), with an overall R² value of 0.93. This result demonstrates the synergistic constraint effect of the dual-channel input on both the deep (gravity anomaly-dominated) and shallow (gravity vertical gradient-dominated) regions, verifying the effectiveness and superiority of this invention.

[0050] Figure 4 This is a graph showing the changes in the loss function values ​​of the training set, validation set, and test set with the number of iterations (Epochs) provided in this embodiment of the invention. The horizontal axis represents the number of training epochs, and the vertical axis represents the loss function value. Figure 4 As can be seen, with the increase in the number of training epochs, the training loss decreases rapidly from an initial high value, then gradually flattens out after about 50 epochs, eventually converging to a lower value. This indicates that the network effectively learns the mapping relationship from the dual-channel input data to the density model on the training set. The validation loss decreases synchronously with the training loss, remaining close to it throughout the training process without significant divergence or an upward trend, indicating that the network has not overfitted. The validation loss fluctuates slightly after convergence, which is a normal phenomenon in deep learning training. After training, the model is evaluated using the test set. The test loss is marked as dots in the graph, and its value is basically consistent with the convergence values ​​of the training and validation losses, further verifying the model's generalization ability.

[0051] The loss curve example shows that the MSE+MAE weighted combined loss function used in this invention can effectively guide network training. The model has achieved good convergence performance on the training set, validation set and test set, and no obvious overfitting or underfitting phenomenon is observed.

[0052] Figure 5This is a visualization example of the joint inversion results. The image shows the inversion results of a random walk test sample, using a 2×2 subplot layout: the top left subplot is the actual 3D density model (3D voxel rendering), with red voxels representing density volumes (ρ=1000kg / m³), and a fixed color scale of 0~1.0g / cm³; the top right subplot is the algorithm-predicted 3D density model (3D voxel rendering), using the same color scale as the actual model for a direct comparison of morphological recovery accuracy. The bottom left subplot is the contour map (unit: mGal) of the actual gravity anomaly Δg, obtained through forward modeling using the kernel matrix G based on the actual density model; the bottom right subplot is the contour map of the gravity anomaly Δg derived from the predicted density model, compared with the bottom left subplot to verify the consistency of the forward modeling. As can be seen from the figure, although the bottom row only shows the gravity anomaly contour lines, the accuracy of the prediction model (top right) stems from the dual-channel input—the Δg channel provides deep constraints (ln term far-field contribution), and the VGG channel provides shallow constraints (arctan term near-field sensitivity). Together, they ensure accurate recovery of both deep and shallow density volumes. The contour lines of the predicted anomaly in the bottom right corner closely match the actual anomaly in the bottom left corner, indicating that the inversion model not only has correct geometric structure but also physically satisfies the observational constraints—something traditional methods cannot simultaneously guarantee.

[0053] In another exemplary embodiment, based on the same inventive concept as the method embodiment, a joint inversion system for gravity anomalies and vertical gradients is provided, comprising: The kernel matrix generation module is used to discretize the underground space into multiple cubic units and calculate the gravity kernel matrix and gravity vertical gradient kernel matrix respectively. The forward modeling module is used to generate density models. It performs forward modeling calculations using the gravity kernel matrix and the gravity vertical gradient kernel matrix to obtain gravity anomaly data and gravity vertical gradient data corresponding to the same density model. The dual-channel input UNet++ network building module is used to build a dual-channel input UNet++ network, which stacks gravity anomaly data and gravity vertical gradient data along the channel dimension into a dual-channel input tensor; A dual-channel input UNet++ network training module is used to train the dual-channel input UNet++ network using a joint inversion loss function; The inversion module is used to input the gravity anomaly data and gravity vertical gradient data to be inverted into the trained dual-channel input UNet++ dense skip connection network to obtain the prediction density model.

[0054] In another exemplary embodiment, based on the same inventive concept as the method embodiment, an electronic device is provided, including a memory and a processor. The memory stores computer instructions that can be executed on the processor. When the processor executes the computer instructions, it performs a gravity anomaly and vertical gradient joint inversion method provided by the embodiment of the present invention.

[0055] The processor may be a single-core or multi-core central processing unit or a specific integrated circuit, or one or more integrated circuits configured to implement the present invention.

[0056] The embodiments of the subject matter and functional operation described in this specification can be implemented in: tangibly embodied computer software or firmware, computer hardware including the structures disclosed in this specification and their structural equivalents, or combinations thereof. Embodiments of the subject matter described in this specification can be implemented as one or more computer programs, i.e., one or more modules of computer program instructions encoded on a tangible, non-transitory program carrier for execution by a data processing device or for controlling the operation of a data processing device. Alternatively or additionally, the program instructions may be encoded on artificially generated propagation signals, such as machine-generated electrical, optical, or electromagnetic signals, which are generated to encode information and transmit it to a suitable receiving device for execution by the data processing device.

[0057] The processing and logic flow described in this specification can be executed by one or more programmable computers that execute one or more computer programs to perform corresponding functions by operating on input data and generating output. The processing and logic flow can also be executed by dedicated logic circuitry—such as FPGAs (Field-Programmable Gate Arrays) or ASICs (Application-Specific Integrated Circuits), and the device can also be implemented as dedicated logic circuitry.

[0058] Suitable processors for executing computer programs include, for example, general-purpose and / or special-purpose microprocessors, or any other type of central processing unit. Typically, the central processing unit receives instructions and data from read-only memory and / or random access memory. The basic components of a computer include a central processing unit for implementing or executing instructions and one or more memory devices for storing instructions and data. Typically, a computer will also include one or more mass storage devices for storing data, such as disks, magneto-optical disks, or optical disks, or the computer will be operatively coupled to such mass storage devices to receive data from or transfer data to them, or both. However, a computer is not required to have such devices. Furthermore, a computer can be embedded in another device, such as a mobile phone, a personal digital assistant (PDA), a mobile audio or video player, a game console, a global positioning system (GPS) receiver, or a portable storage device such as a universal serial bus (USB) flash drive, to name a few.

[0059] It should be understood that each block in a flowchart or block diagram can represent a module, segment, or portion of code, which contains one or more executable instructions for implementing the specified logical function. It should also be noted that in some alternative implementations, the functions marked in the blocks may occur in a different order than those shown in the figures. For example, two consecutive blocks may actually be executed substantially in parallel, or they may sometimes be executed in reverse order, depending on the functions involved. It should also be noted that each block in a block diagram and / or flowchart, and combinations of blocks in block diagrams and / or flowcharts, can be implemented using a dedicated hardware-based system that performs the specified function or action, or using a combination of dedicated hardware and computer instructions.

[0060] The above detailed embodiments are a description of the present invention. It should not be considered that the specific embodiments of the present invention are limited to these descriptions. For those skilled in the art, several simple deductions and substitutions can be made without departing from the concept of the present invention, and all of these should be considered to fall within the protection scope of the present invention.

Claims

1. A method for joint inversion of gravity anomalies and vertical gradients, characterized in that, Includes the following steps: S1. Discretize the underground space into multiple cubic elements and calculate the gravity kernel matrix and gravity vertical gradient kernel matrix respectively; S2. Generate a density model, and perform forward modeling using the gravity kernel matrix and the gravity vertical gradient kernel matrix respectively to obtain gravity anomaly data and gravity vertical gradient data corresponding to the same density model; S3. Construct a dual-channel input UNet++ network to stack gravity anomaly data and gravity vertical gradient data along the channel dimension into a dual-channel input tensor; S4. The dual-channel input UNet++ network is trained using a joint inversion loss function; S5. Input the gravity anomaly data and gravity vertical gradient data to be inverted into the trained dual-channel input UNet++ dense skip connection network to obtain the prediction density model.

2. The method for joint inversion of gravity anomalies and vertical gradients according to claim 1, characterized in that, The dual-channel input UNet++ network adopts a dense skip connection structure, in which the encoder downsamples layer by layer to extract multi-scale features, and the decoder fuses the features of each layer through dense skip connections. The encoder has two input channels in its first convolutional layer, which are used to simultaneously receive gravity anomaly data and gravity vertical gradient data.

3. The method for joint inversion of gravity anomalies and vertical gradients according to claim 1, characterized in that, The gravity kernel matrix is ​​calculated based on the Nagy rectangular gravity anomaly analytical formula, and the kernel function contains ln terms and arctan terms; the gravity vertical gradient kernel matrix is ​​obtained by taking the z-direction partial derivative of the gravity formula, and the kernel function contains only arctan terms.

4. The method for joint inversion of gravity anomalies and vertical gradients according to claim 2, characterized in that, The dual-channel input tensor has a shape of (2, H, W), with one channel inputting gravity anomaly data and the other channel inputting gravity vertical gradient data; the encoder's first layer convolutional kernel size is (2, 3, 3), and the kernel weights are automatically learned through backpropagation to obtain the optimal fusion weights for gravity anomaly data and gravity vertical gradient data.

5. The method for joint inversion of gravity anomalies and vertical gradients according to claim 1, characterized in that, The joint inversion loss function is a weighted combination of mean squared error (MSE) and mean absolute error (MAE): L = α·MSE + β·MAE, where α and β are weighting coefficients, the mean square error (MSE) is used to eliminate large deviations to restore the global shape, and the mean absolute error (MAE) is used to reduce small deviations to optimize local details.

6. The method for joint inversion of gravity anomalies and vertical gradients according to claim 1, characterized in that, During training in step S4, the z-score of the dual-channel input is independently normalized and the normalization parameters of each channel are persistently saved.

7. The method for joint inversion of gravity anomalies and vertical gradients according to claim 1, characterized in that, Gravity anomaly data and gravity vertical gradient data corresponding to the same density model are paired and assembled according to their numbers to generate joint data, ensuring that the two types of observation data share the same source model.

8. The method for joint inversion of gravity anomalies and vertical gradients according to claim 1, characterized in that, After obtaining the predicted density model, the gravity kernel matrix and the gravity vertical gradient kernel matrix are loaded. Based on the predicted density model, the predicted gravity anomaly and the predicted gravity vertical gradient are calculated in forward modeling. They are then compared and verified with the observation data to achieve joint evaluation of the dual kernel matrix.

9. A joint inversion system for gravity anomalies and vertical gradients, characterized in that, include: The kernel matrix generation module is used to discretize the underground space into multiple cubic units and calculate the gravity kernel matrix and gravity vertical gradient kernel matrix respectively. The forward modeling module is used to generate density models. It performs forward modeling calculations using the gravity kernel matrix and the gravity vertical gradient kernel matrix to obtain gravity anomaly data and gravity vertical gradient data corresponding to the same density model. The dual-channel input UNet++ network building module is used to build a dual-channel input UNet++ network, which stacks gravity anomaly data and gravity vertical gradient data along the channel dimension into a dual-channel input tensor; A dual-channel input UNet++ network training module is used to train the dual-channel input UNet++ network using a joint inversion loss function; The inversion module is used to input the gravity anomaly data and gravity vertical gradient data to be inverted into the trained dual-channel input UNet++ dense skip connection network to obtain the prediction density model.

10. An electronic device comprising a memory and a processor, wherein the memory stores computer instructions executable by the processor, characterized in that, When the processor executes computer instructions, it performs a joint inversion method for gravity anomalies and vertical gradients as described in any one of claims 1-8.

Images