Real-time regional gravity field modeling method and system based on deep neural network

Abstract

The invention relates to the technical field of deep learning, and discloses a real-time regional gravity field modeling method and system based on a deep neural network, and the method comprises the following steps: obtaining ground gravity measurement data, and calculating the original ground gravity anomaly of ground discrete points; a global gravity field model is introduced, and long-wave gravity abnormal components of ground discrete points are obtained in real time; a deep full-connection neural network FC-DNN is introduced; carrying out real-time calculation to obtain short-wave gravity abnormal components of ground discrete points; introducing a convolutional neural network (CNN); performing real-time prediction to obtain residual gravity abnormal components of ground discrete points; and fusing the residual gravity abnormal component, the short-wave gravity abnormal component and the long-wave gravity abnormal component of the ground discrete points to complete modeling of the regional ground gravity field model. The method has the advantages of being high in training efficiency and reasoning speed, and the problem that in existing gravity field modeling, the calculation efficiency is insufficient is solved.

Description

Technical Field

[0001] The present invention relates to the field of deep learning technology, and more specifically, to a real-time regional gravity field modeling method and system based on deep neural network. Background Art

[0002] A high-precision local gravity field is the foundation and basis for regional geoid refinement, engineering surveying, navigation and positioning, geophysical exploration, and geological structure analysis. Therefore, research on how to establish more accurate and reliable local gravity field models to serve the needs of national economic development and scientific research is of great significance. Refining multi-source data preprocessing methods and improving gravity field modeling technology have been hot topics in recent years in the field of gravity field modeling.

[0003] Modern local gravity field modeling is primarily based on the removal-recovery theory, combining high-precision ground gravity data, digital terrain models, and GPS leveling. The specific process includes: 1) preprocessing raw gravity anomaly data acquired by gravimeters or satellites, including correction and filtering; 2) removing the long-wavelength gravity effect and the topographic gravity effect from the observed gravity anomaly; 3) mathematically refining the residual gravity anomaly; and 4) restoring the long-wavelength gravity effect and the topographic gravity effect based on the refined residual gravity anomaly to complete the local gravity field modeling.

[0004] In the removal-restore framework, terrain gravity effects are typically calculated using high-precision terrain data and the law of universal gravitation. However, currently used calculation methods each have limitations. Spectral domain methods suffer from series divergence, while spatial domain methods suffer from low computational efficiency. Among the most widely used methods for refining residual gravity anomalies, the radial basis function method and the least squares method are both involved in fusing multi-source data, leading to difficulties in data weighting and large computational load. The massive computational effort required in these two steps results in inefficient and poorly maneuverable regional gravity field modeling, hindering the further development of gravity field models. For example, in space science, the motion of near-Earth objects is affected by gravity, and Earth's gravity field parameters are important mechanical parameters in orbital dynamics models. Rapid, real-time calculation of gravity field parameters would greatly enhance the utility of gravity field models.

[0005] In view of the computational efficiency problems existing in current gravity field modeling, how to invent a real-time regional gravity field modeling method with flexible input and powerful nonlinear fitting capabilities is a technical problem that urgently needs to be solved in this technical field. Summary of the Invention

[0006] In order to solve the problem of insufficient computational efficiency in gravity field modeling, the present invention provides a real-time regional gravity field modeling method and system based on deep neural network, which has the characteristics of high training efficiency and fast inference speed.

[0007] In order to achieve the above-mentioned purpose of the present invention, the technical solutions adopted are as follows:

[0008] A real-time regional gravity field modeling method based on a deep neural network includes the following specific steps:

[0009] Obtain ground gravity measurement data and calculate the original ground gravity anomaly of its discrete points;

[0010] The global gravity field model is introduced to calculate the long-wave gravity anomaly components of discrete points on the ground in real time;

[0011] A terrain data sample set is obtained; a deep fully connected neural network (FC-DNN) is introduced; the sample locations and terrain features in the terrain data sample set are used as input, and the shortwave gravity anomaly is used as the supervision target to train the FC-DNN to learn the mapping relationship between terrain features and shortwave gravity anomaly. The trained FC-DNN is used as a fast shortwave gravity field solution model, and the shortwave gravity anomaly component is obtained in real time based on the terrain features of discrete ground points.

[0012] The long-wave and short-wave gravity anomaly components of discrete ground points are removed from their original ground gravity anomalies to obtain their residual gravity anomalies; a multi-source data sample set is constructed, and a CNN is trained to model the nonlinear relationship between the multi-source data and the residual gravity anomaly; the trained CNN is used as a fast solution model for the residual gravity anomaly, and the residual gravity anomaly components are obtained by real-time prediction based on the multi-source data of discrete ground points;

[0013] The residual gravity anomaly component of discrete ground points, the short-wave gravity anomaly component, and the long-wave gravity anomaly component are integrated to complete the modeling of the regional ground gravity field model.

[0014] Preferably, the original ground gravity anomaly is calculated by the following specific steps:

[0015] Obtain ground gravity measurement data for each sample, including longitude, latitude, surface elevation, height of discrete ground points, and the difference between the height of discrete ground points and the surface elevation;

[0016] Perform preprocessing to eliminate systematic errors and outliers in ground gravity measurement data;

[0017] The normal gravity is subtracted from the preprocessed gravity measurement data to obtain the original ground gravity anomaly of its discrete points.

[0018] Furthermore, the global gravity field model is introduced to construct the long-wave gravity anomaly components of discrete points on the ground. The specific steps are as follows:

[0019] Selecting the global gravity field model, the gravitational potential V at a point P outside the Earth is defined by its geocentric distance r, geocentric colatitude θ and longitude λ:

[0020]

[0021] Where GM is the Earth's gravitational constant, a is the fully normalized spherical harmonic coefficient Related scale factors, in addition:

[0022]

[0023] in, is a fully normalized Legendre function of the first kind, with order n and degree |m|. The first-order term (n=1) in the spherical harmonic expansion represents the offset of the center of gravity.

[0024] The disturbance potential T at point P(r,θ,λ) is defined as the difference between the actual gravity potential of the Earth and the normal gravity potential;

[0025] The perturbation potential is expressed by spherical harmonic expansion:

[0026]

[0027] The zero-order term in the formula reflects the gravitational effect produced by the total mass of the Earth. If the mass of the reference ellipsoid is equal to the actual mass of the Earth, this term is zero.

[0028] The surface free-air gravity anomaly Δg is defined as the actual gravity acceleration at point P minus the normal gravity acceleration at the corresponding geoid point Q:

[0029] Δg=|g P |-|γ Q |

[0030] Assuming that the gravitational potential on the geoid is equal to the normal gravitational potential on the reference ellipsoid, the basic relationship between gravity anomaly and disturbance potential is obtained:

[0031]

[0032] in represents the radial derivative of the perturbation potential, is a correction term related to the distance, and the three ε terms represent three types of error corrections;

[0033] Let Δg c Represents the gravity anomaly after all systematic corrections have been applied:

[0034]

[0035] get:

[0036]

[0037] The long-wave gravity anomaly component is further obtained:

[0038]

[0039] in

[0040]

[0041] From this we get The spherical harmonic expansion expression of , whose infinite series in the external space of the Brillouin sphere has convergence.

[0042] Furthermore, we can obtain a sample set of terrain data by the following steps:

[0043] Define source quality models; collect high-resolution terrain feature data and extract terrain features including surface elevation data;

[0044] The high-resolution terrain feature data is filtered to obtain reference terrain data, and the difference between the reference terrain data and the high-resolution terrain feature data is used as the residual terrain model RTM;

[0045] Randomly select discrete points from the reference terrain data, record the sample position of each discrete point including longitude, latitude and altitude, and calculate its shortwave gravity anomaly;

[0046] The RTM mass distribution is divided into several regions according to the distance from the calculation point, and the shortwave gravity anomaly is approximately calculated for each region using geometric bodies;

[0047] The terrain data sample set is obtained by fusing terrain features, shortwave gravity anomalies, and sample locations.

[0048] Furthermore, the FC-DNN includes an input layer, L hidden layers, and an output layer. The FC-DNN is trained in an end-to-end manner to learn the mapping relationship between terrain features and shortwave gravity anomalies, specifically:

[0049] Assume that there are N samples of terrain data and adopt batch training method. Set the number of samples in each batch to k. According to N / k=q, divide the data into q batches. Each batch performs forward operation to obtain the loss function, and then backpropagates to update the weights. This process is called one iteration.

[0050] The forward propagation of each layer is defined as:

[0051] H (0)=X

[0052] Z (l) =W (l) H (l-1) +B (l)

[0053] H (l) =σ(Z (l) )

[0054] Among them, X is the matrix of input samples, with dimension [m×k], W (l) is the weight, B (l) is the bias, which is broadcast to k columns during calculation to match the batch size;

[0055] The batch loss function for training is:

[0056]

[0057] Where Y is the true label matrix with dimension [n L ×k],H (L) Is the output matrix of the output layer, with dimension [n L ×k],||·|| F represents the Frobenius norm;

[0058] During the backpropagation process, the error term of the output layer is calculated as:

[0059]

[0060] Where ⊙ represents the element-by-element multiplication of the matrix;

[0061] For the hidden layer (l = L-1, L-2, ..., 1), the recursive formula of the error term is:

[0062] Δ (l) =(W (l+1) ) T Δ (l+1) ⊙σ′(Z (l) )

[0063] The batch gradient is calculated as:

[0064]

[0065] Among them 1 k is a k-dimensional all-one column vector used to sum the gradients over the batch dimension;

[0066] The parameter update formula is:

[0067]

[0068] In batch training, all batches undergo one iteration and complete one cycle, and training continues until the set threshold is reached to obtain the trained FC-DNN.

[0069] Furthermore, the FC-DNN training also introduces a terrain regularization method to further improve the generalization ability and prediction accuracy of the network; the loss function of the regularization method is:

[0070]

[0071] Where a is a trainable array with a dimension of [1×k]. The mean and variance used to initialize it are the mean and variance of the gravity anomaly and terrain ratio in the training area. Predict gravity anomalies for the training area, h is H (L) The corresponding surface elevation value.

[0072] Furthermore, the structure of CNN includes:

[0073] Feature extraction part: includes a convolution layer, batch normalization layer, activation function layer, and pooling layer;

[0074] Position information fusion part: flatten the feature map output by the feature extraction part and then splice it with the position information of the discrete points on the ground;

[0075] Regression prediction part: includes at least one fully connected layer, which is used to map the fused feature vector into the residual gravity anomaly prediction value.

[0076] Furthermore, the training process of CNN is as follows:

[0077] Assume that the input of the lth convolutional layer in the feature extraction part is h (l-1) , the convolution kernel is K (l) , bias is b (l) , define the output of the convolutional layer as:

[0078] h (l) =f(Conv(h (l-1) ,K (l) )+b(l))

[0079] Where f(·) is the activation function and Conv(·) is the convolution operation, which is defined as:

[0080]

[0081] Among them O i,j is the value of the output feature map at position (i, j), k is the convolution kernel size, and b is the bias term;

[0082] The batch normalization layer is defined as:

[0083]

[0084] Among them, μ B and is the mean and variance of the batch, γ and β are learnable scaling and translation parameters, is a constant;

[0085] The output of the batch normalization layer passes through the Tanh activation function layer and then the pooling layer to further reduce the feature map size by downsampling:

[0086]

[0087] Where R is the pooling window, and features are extracted layer by layer by stacking multiple convolutional layers and pooling layers;

[0088] In the position information fusion part, the extracted feature map is flattened and then spliced ​​with the position of the discrete points on the ground to form a complete feature vector;

[0089] The complete feature vector is flattened into a one-dimensional vector, and the fully connected layer is connected to the output layer for regression task;

[0090] Introducing the weighted loss function:

[0091] Construct a gradient descent method to minimize the loss function E:

[0092]

[0093] Among them, w i is the error weight, according to y i The distribution definition of y i When the distribution is uneven, values ​​in intervals with fewer numerical distributions are assigned greater weights;

[0094] Reshape the error into the shape of the feature map:

[0095]

[0096] For the max pooling layer, the gradient is only passed to the location of the maximum value in the pooling window:

[0097]

[0098] Where |R| is the size of the pooling window;

[0099] The gradient of the batch normalization layer is calculated as:

[0100]

[0101] The gradient of the loss function to the convolution kernel in the output gradient of the convolution layer l is:

[0102]

[0103] The gradient with respect to the bias is:

[0104]

[0105] The gradient of the input feature map is:

[0106]

[0107] Calculated As the error input of the previous layer, continue backpropagation training;

[0108] Update the network parameters using gradient descent:

[0109]

[0110] The input multi-source data sample set is divided into a training set and a test set. The training set is divided into multiple batches and trained for several cycles. Some samples are selected from the training set as a validation set. The network accuracy is verified after each cycle of training. At the same time, the hyperparameters are adjusted to ensure that the accuracy meets expectations and there is no overfitting. The trained CNN is obtained.

[0111] Furthermore, to address the gradient vanishing and network degradation problems faced by CNN, at least one residual block is introduced in the feature extraction part, and its structure is as follows:

[0112] y=σ(F(x,{W i})+x)

[0113] Where σ represents the nonlinear transformation composed of the convolution layer, batch normalization layer and activation function; x is the feature map extracted by the feature extraction part, and y is the complete feature map output by the residual block;

[0114] For the residual block, the error backpropagation is expressed as:

[0115]

[0116] h (l) =f(Conv(h (l-1) ,K (l) )+b(l)).

[0117] A real-time regional gravity field modeling system based on deep neural network, including gravity measurement module, long-wave gravity anomaly extraction module, short-wave gravity anomaly extraction module, residual gravity anomaly extraction module, and fusion reconstruction module;

[0118] The gravity measurement module is used to obtain ground gravity measurement data and calculate the original ground gravity anomaly of its discrete points on the ground;

[0119] The long-wave gravity anomaly extraction module is used to introduce the global gravity field model and calculate the long-wave gravity anomaly components of discrete points on the ground in real time;

[0120] The shortwave gravity anomaly extraction module is used to obtain a terrain data sample set; a deep fully connected neural network (FC-DNN) is introduced; the sample locations and terrain features in the terrain data sample set are used as input, and the shortwave gravity anomaly therein is used as a supervision target to train the FC-DNN to learn the mapping relationship between the terrain features and the shortwave gravity anomaly; the trained FC-DNN is used as a shortwave gravity field fast solution model, and the shortwave gravity anomaly components are obtained based on the terrain features of discrete points on the ground in real time.

[0121] The residual gravity anomaly extraction module is used to remove the long-wave and short-wave gravity anomaly components of the ground discrete points from their original ground gravity anomalies to obtain their residual gravity anomalies; construct a multi-source data sample set, train a CNN to model the nonlinear relationship between the multi-source data and the residual gravity anomaly; use the trained CNN as a rapid solution model for the residual gravity anomaly, and obtain its residual gravity anomaly components based on real-time prediction of the multi-source data of the ground discrete points;

[0122] The fusion reconstruction module is used to fuse the residual gravity anomaly component, the short-wave gravity anomaly component, and the long-wave gravity anomaly component of the ground discrete points to complete the modeling of the regional ground gravity field model.

[0123] The beneficial effects of the present invention are as follows:

[0124] The present invention discloses a real-time regional gravity field modeling method based on deep neural networks. In response to the computational efficiency issues existing in current gravity field modeling, a modeling method combining deep fully connected neural networks and convolutional neural networks is proposed. Through deep learning technology combined with the removal-recovery method, effective fusion and rapid calculation of multi-source data are achieved. The method can be used for the real-time construction of regional gravity field models, providing high-precision and high-efficiency gravity field data support for fields such as geodesy, geophysical exploration, and space science. BRIEF DESCRIPTION OF THE DRAWINGS

[0125] Figure 1 It is a flow chart of a real-time regional gravity field modeling method based on a deep neural network of the present invention.

[0126] Figure 2 This is a specific flow chart of a real-time regional gravity field modeling method based on a deep neural network in Example 2 of the present invention.

[0127] Figure 3 This is a schematic diagram of the training process of FC-DNN, a real-time regional gravity field modeling method based on deep neural network in Example 2 of the present invention.

[0128] Figure 4 This is a schematic diagram of the CNN training process of a real-time regional gravity field modeling method based on a deep neural network in Example 2 of the present invention. DETAILED DESCRIPTION

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

[0130] Example 1

[0131] like Figure 1 As shown, a real-time regional gravity field modeling method based on deep neural network includes the following specific steps:

[0132] Obtain ground gravity measurement data and calculate the original ground gravity anomaly of its discrete points;

[0133] The global gravity field model is introduced to calculate the long-wave gravity anomaly components of discrete points on the ground in real time;

[0134] A terrain data sample set is obtained; a deep fully connected neural network (FC-DNN) is introduced; the sample locations and terrain features in the terrain data sample set are used as input, and the shortwave gravity anomaly is used as the supervision target to train the FC-DNN to learn the mapping relationship between terrain features and shortwave gravity anomaly. The trained FC-DNN is used as a fast shortwave gravity field solution model, and the shortwave gravity anomaly component is obtained in real time based on the terrain features of discrete ground points.

[0135] The long-wave and short-wave gravity anomaly components of discrete ground points are removed from their original ground gravity anomalies to obtain their residual gravity anomalies; a multi-source data sample set is constructed, and a CNN is trained to model the nonlinear relationship between the multi-source data and the residual gravity anomaly; the trained CNN is used as a fast solution model for the residual gravity anomaly, and the residual gravity anomaly components are obtained by real-time prediction based on the multi-source data of discrete ground points;

[0136] The residual gravity anomaly component of discrete ground points, the short-wave gravity anomaly component, and the long-wave gravity anomaly component are integrated to complete the modeling of the regional ground gravity field model.

[0137] Example 2

[0138] like Figure 2 As shown, the real-time regional gravity field modeling method based on deep neural network in this embodiment has the following specific steps:

[0139] Obtain ground gravity measurement data within the regional range, set the location information of discrete points on the ground, and perform preprocessing to eliminate systematic errors and outliers, and calculate the original ground gravity anomaly;

[0140] Select an appropriate global gravity field model and construct the long-wave gravity anomaly components of discrete ground points based on spherical harmonic analysis and spherical harmonic synthesis methods;

[0141] By integrating high-resolution terrain feature data, using the residual terrain model (RTM) and forward modeling technology to solve shortwave gravity anomalies and construct a shortwave gravity anomaly sample set, a deep fully connected neural network (FC-DNN) is introduced. Using sample location information and terrain elevation as input and gravity anomalies caused by high-frequency terrain as supervision targets, the FC-DNN is trained to learn the mapping relationship between terrain features and shortwave gravity anomalies, replacing the classical integral model to achieve efficient prediction of local shortwave gravity anomalies. On this basis, a terrain regularization method is introduced to further improve the network's generalization ability and prediction accuracy.

[0142] The long-wave and short-wave gravity anomalies of the obtained ground discrete points are removed from the ground gravity anomaly to obtain the residual gravity anomaly;

[0143] A dataset was constructed using multi-source data, such as landforms and crustal thickness, and the nonlinear relationship between this data and ground gravity anomalies was modeled using a convolutional neural network (CNN). The CNN extracts key patterns from this heterogeneous multi-source information through spatial slicing, scale unification, and feature fusion, improving the model's ability to predict local gravity anomalies. Furthermore, a weighted loss function was introduced to adjust for the uneven distribution of residual gravity anomalies, further improving the balance and accuracy of the network's predictions.

[0144] The residual gravity anomaly component of discrete ground points predicted by the trained CNN, the short-wave gravity anomaly calculated by the trained FC-DNN, and the long-wave gravity anomaly provided by the spherical harmonic model are integrated to complete the reconstruction of the regional ground gravity field model, achieving the goal of high-precision and high-efficiency refinement of the local gravity field.

[0145] In a specific embodiment, the original ground gravity anomaly is calculated by the following steps:

[0146] Obtain ground gravity measurement data for each sample, including longitude, latitude, surface elevation, height of discrete ground points, and the difference between the height of discrete ground points and the surface elevation;

[0147] Perform preprocessing to eliminate systematic errors and outliers in ground gravity measurement data;

[0148] The normal gravity is subtracted from the preprocessed gravity measurement data to obtain the original ground gravity anomaly of its discrete points.

[0149] In a specific embodiment, a global gravity field model is introduced to construct long-wave gravity anomaly components of discrete points on the ground. The specific steps are as follows:

[0150] Selecting the global gravity field model, the gravitational potential V at a point P outside the Earth is defined by its geocentric distance r, geocentric colatitude θ and longitude λ:

[0151]

[0152] Where GM is the Earth's gravitational constant, a is the fully normalized spherical harmonic coefficient The relevant scale factor is usually taken as the equatorial radius of the mean Earth ellipsoid, and in addition:

[0153]

[0154] in, is a fully normalized Legendre function of the first kind, with order n and degree |m|. The first-order term (n=1) in the spherical harmonic expansion represents the offset of the center of gravity. In the geocentric coordinate system, the geometric parameters of the reference ellipsoid already reflect the overall shape of the Earth, and the offset of the center of gravity in the coordinate system is zero, so the first-order term is ignored.

[0155] The perturbation potential T at a point P(r,θ,λ) is defined as the difference between the actual gravity potential of the Earth and the normal gravity potential; where the normal gravity potential refers to the gravity potential of the reference ellipsoid at that point. In this case, the perturbation potential is expressed by spherical harmonic expansion:

[0156]

[0157] The zero-order term in the formula reflects the gravitational effect produced by the total mass of the Earth. If the mass of the reference ellipsoid is equal to the actual mass of the Earth, this term is zero.

[0158] The surface free-air gravity anomaly Δg is defined as the actual gravity acceleration at point P minus the normal gravity acceleration at the corresponding geoid point Q:

[0159] Δg=|g P |-|γ Q |

[0160] Assuming that the gravitational potential on the geoid is equal to the normal gravitational potential on the reference ellipsoid, the basic relationship between gravity anomaly and disturbance potential is obtained:

[0161]

[0162] in represents the radial derivative of the perturbation potential, is a correction term related to the distance, and the three ε terms represent three types of error corrections;

[0163] Let Δg c Represents the gravity anomaly after all systematic corrections have been applied:

[0164]

[0165] get:

[0166]

[0167] The long-wave gravity anomaly component is further obtained:

[0168]

[0169] Gravity anomaly Δg c (r,θ,λ) itself is not a harmonic function, but The space outside the mass is harmonious, where

[0170]

[0171] From this we get The spherical harmonic expansion expression of , whose infinite series in the external space of the Brillouin sphere has convergence.

[0172] In a specific embodiment, Figure 3 As shown, the terrain data sample set is obtained. The specific steps are:

[0173] The normal density is 2.67g / cm 3 The terrain quality is defined as the source quality model; MERIT DEM (Multi-Error-Removed Improved-Terrain DEM) is used as high-resolution terrain feature data, and terrain features including surface elevation data are extracted;

[0174] The spatial resolution of MERIT DEM is 3″. The high-resolution terrain feature data are filtered using a filter with a window size of 100 × 100 to obtain the reference terrain data. The difference between the reference terrain data and the high-resolution terrain feature data is used as the residual terrain model (RTM);

[0175] Randomly select discrete points from the reference terrain data and record the location information of each discrete point, including longitude, latitude, surface elevation, height of the ground discrete point, and the difference between the height of the ground discrete point and the surface elevation; the terrain quality is called the source quality model, and the gravitational potential V generated by the random discrete point at point P is p Expressed as:

[0176]

[0177] Among them, G is the gravitational constant, M and m i are the mass of the source mass model and the mass of a single integration element, dm and dm i M and m respectively i The mass element, r p and m i Calculate points P to M and m respectively i Euclidean distance of the geometric center;

[0178] The coordinate axis is established with the calculation point P as the origin, with the x axis pointing to the north of the earth, the y axis pointing to the east of the earth, and the z axis pointing to the zenith, and the shortwave gravity anomaly Δg at point P is obtained. z :

[0179]

[0180] Thus, the shortwave gravity anomalies of all discrete points are obtained;

[0181] The RTM mass distribution is divided into three regions based on the distance from the calculation point and approximated using geometric bodies such as polyhedrons, prisms, Tesseroids, and point masses. The terrain gravity effect within 0.03° of the calculation point is calculated using a prismatic model with a source mass model spatial resolution of 3″; the range from 0.03° to 0.80° is calculated using a Tesseroid model with a source mass model spatial resolution of 3″; and the range from 0.80° to 2.00° is calculated using a point mass model with a source mass model spatial resolution of 30″. Shortwave gravity anomalies are approximated for each region using geometric bodies.

[0182] The terrain data sample set is obtained by fusing terrain features, shortwave gravity anomalies, and sample locations.

[0183] In a specific embodiment, the FC-DNN includes an input layer, L hidden layers, and an output layer. During training, the FC-DNN learns the mapping relationship between terrain features and shortwave gravity anomalies in an end-to-end manner, specifically:

[0184] In this embodiment, the input layer of the FC-DNN specifically includes 5 neurons and 6 hidden layers; each hidden layer contains 256 neurons, and all hidden layers are fully connected; the output layer contains 1 neuron, corresponding to the only output variable - shortwave gravity anomaly; the activation function of the hidden layer is selected as Tanh, the activation function of the output layer is selected as Linear, and the optimizer is selected as Adam.

[0185] Assume that there are N samples of terrain data and adopt batch training method. Set the number of samples in each batch to k. According to N / k=q, divide the data into q batches. Each batch performs forward operation to obtain the loss function, and then backpropagates to update the weights. This process is called one iteration.

[0186] The forward propagation of each layer is defined as:

[0187] H (0) =X

[0188] Z (l) =W (l) H (l-1) +B (l)

[0189] H (l) =σ(Z (l) )

[0190] Among them, X is the matrix of input samples, with dimension [m×k], W (l) is the weight, B (l) is the bias, which is broadcast to k columns during calculation to match the batch size;

[0191] The batch loss function for training is:

[0192]

[0193] Where Y is the true label matrix with dimension [n L ×k],H (L) Is the output matrix of the output layer, with dimension [n L ×k],||·|| F represents the Frobenius norm;

[0194] During the backpropagation process, the error term of the output layer is calculated as:

[0195]

[0196] Where ⊙ represents the element-by-element multiplication of the matrix;

[0197] For the hidden layer (l = L-1, L-2, ..., 1), the recursive formula of the error term is:

[0198] Δ (l) =(W (l+1) )TΔ (l+1) ⊙σ′(Z (l) )

[0199] The batch gradient is calculated as:

[0200]

[0201]

[0202] Among them 1 k is a k-dimensional all-one column vector used to sum the gradients over the batch dimension;

[0203] The parameter update formula is:

[0204]

[0205] In batch training, all batches undergo one iteration and complete one cycle, and training continues until the set threshold is reached to obtain the trained FC-DNN.

[0206] Furthermore, the FC-DNN training also introduces a terrain regularization method to further improve the generalization ability and prediction accuracy of the network; the loss function of the regularization method is:

[0207]

[0208] Where a is a trainable array with a dimension of [1×k]. The mean and variance used to initialize it are the mean and variance of the gravity anomaly and terrain ratio in the training area. Predict gravity anomalies for the training area, h is H (L) The corresponding surface elevation value.

[0209] In this embodiment, after training the FC-DNN, its accuracy was also tested. 80% of the training samples were selected as the training set, and the learning rate was set to 0.0005, the batch size was 4096, and the training cycle was 3000 for training. 20% of the samples were randomly selected as the test set to test the network performance. The test set does not participate in the parameter optimization process.

[0210] The FC-DNN that passed the accuracy test is used as a fast solution model for the shortwave gravity field to quickly solve the shortwave gravity anomalies of discrete points on the surface.

[0211] In this embodiment, the pseudo code of the FC-DNN training algorithm is:

[0212] Input: training set D = {(x i ,y i )} i＝1,…,N ;

[0213] Batch size k;

[0214] Learning rate η

[0215] The maximum training period max_epochs.

[0216] process:

[0217]

[0218]

[0219] Output: trained neural network parameters W(l) and B (l) .

[0220] Example 3

[0221] In this embodiment, the method for constructing the multi-source data sample set is specifically as follows: collecting multi-source data including but not limited to DEM, reference topography, surface cover, crust thickness, sediment thickness, surface cover type of Moho surface depth, crust thickness, etc., performing pre-processing such as splicing, fusion, and interpolation on the multi-source data, and according to the position of the residual gravity anomaly, each discrete point corresponds to a group of samples of grid data with a diameter of 0.6° around it as a group of samples to obtain a multi-source data sample set.

[0222] In a specific embodiment, Figure 4 As shown, the structure of CNN includes:

[0223] Feature extraction part: includes a convolution layer, a batch normalization layer, an activation function layer, a pooling layer, and two residual blocks;

[0224] In this example, the input features of a multi-source data sample set, including DEM, reference topography, surface cover, crust thickness, sediment thickness, and Moho depth, are extracted through a 3×3 convolutional layer, outputting 16 feature maps. These maps are then nonlinearly transformed using a batch normalization layer and a Tanh activation function, and finally reduced in size using a max pooling layer. The resulting feature maps are then passed through two residual blocks for deep feature extraction. Each residual block contains two 3×3 convolutional layers, paired with a batch normalization layer and a Tanh activation function, with the number of channels expanded to 32 and 64, respectively. A random dropout layer is added during the feature extraction process, with a dropout rate of 0.2.

[0225] In this embodiment, when the number of input and output channels of the residual block does not match, 1×1 convolution is used to perform channel adjustment to ensure that the residual connection is correctly implemented;

[0226] Position information fusion part: The feature map obtained by the two residual blocks is flattened and then spliced ​​with the position code including longitude, latitude, and elevation to form a complete feature vector;

[0227] The regression prediction part consists of a three-layer fully connected network. The first layer reduces the concatenated feature vector to 256 dimensions, and the second layer further reduces it to 128 dimensions. The final output is a single predicted value - the residual gravity anomaly.

[0228] In this embodiment, the CNN training process is specifically as follows:

[0229] Assume that the input of the lth convolutional layer in the feature extraction part is h (l-1) , the convolution kernel is K (l) , bias is b(l) , define the output of the convolutional layer as:

[0230] h (l) =f(Conv(h (l-1) ,K (l) )+b(l))

[0231] Where f(·) is the activation function and Conv(·) is the convolution operation, which is defined as:

[0232]

[0233] Among them O i,j is the value of the output feature map at position (i, j), k is the convolution kernel size, and b is the bias term;

[0234] The batch normalization layer is defined as:

[0235]

[0236] Among them, μ B and is the mean and variance of the batch, γ and β are learnable scaling and translation parameters, is a constant;

[0237] The output of the batch normalization layer passes through the Tanh activation function layer and then the pooling layer to further reduce the feature map size by downsampling:

[0238]

[0239] Where R is the pooling window, and features are extracted layer by layer by stacking multiple convolutional layers and pooling layers;

[0240] In the position information fusion part, the extracted feature map is flattened and then spliced ​​with the position of the discrete points on the ground to form a complete feature vector;

[0241] The complete feature vector is flattened into a one-dimensional vector, and the fully connected layer is connected to the output layer for regression task;

[0242] Introducing the weighted loss function:

[0243] Construct a gradient descent method to minimize the loss function E:

[0244]

[0245] Among them, w i is the error weight, according to y i The distribution definition of y i When the distribution is uneven, values ​​in intervals with fewer numerical distributions are assigned greater weights;

[0246] Reshape the error into the shape of the feature map:

[0247]

[0248] For the max pooling layer, the gradient is only passed to the location of the maximum value in the pooling window:

[0249]

[0250] Where |R| is the size of the pooling window;

[0251] The gradient of the batch normalization layer is calculated as:

[0252]

[0253] The gradient of the loss function to the convolution kernel in the output gradient of the convolution layer l is:

[0254]

[0255] The gradient with respect to the bias is:

[0256]

[0257] The gradient of the input feature map is:

[0258]

[0259] Calculated As the error input of the previous layer, continue backpropagation training;

[0260] Update the network parameters using gradient descent:

[0261]

[0262] The input multi-source data sample set is divided into a training set and a test set. The training set is divided into multiple batches and trained for several cycles. Some samples are selected from the training set as a validation set. The network accuracy is verified after each cycle of training. At the same time, the hyperparameters are adjusted to ensure that the accuracy meets expectations and there is no overfitting. The trained CNN is obtained.

[0263] In this embodiment, hyperparameters are selected to train the network. Preferably, the input data is split into a training set and a test set at an 8:2 ratio. The training set is divided into multiple batches (Batch Size), and multiple epochs of training are performed. A portion of the training set samples are selected as a validation set. After each epoch of training, the network accuracy is verified, and the hyperparameters are adjusted to ensure that the accuracy meets expectations and there is no overfitting. Specifically, the batch size is set to 128 and the epoch is set to 30.

[0264] In a specific embodiment, to address the gradient vanishing and network degradation problems faced by CNN, at least one residual block is introduced in the feature extraction part, and its structure is as follows:

[0265] y=σ(F(x,{W i})+x)

[0266] Where σ represents the nonlinear transformation composed of the convolution layer, batch normalization layer and activation function; x is the feature map extracted by the feature extraction part, and y is the complete feature map output by the residual block;

[0267] For the residual block, the error backpropagation is expressed as:

[0268]

[0269] h (l) =f(Conv(h (l-1) ,K (l) )+b(l)).

[0270] Example 4

[0271] A real-time regional gravity field modeling system based on deep neural network, including gravity measurement module, long-wave gravity anomaly extraction module, short-wave gravity anomaly extraction module, residual gravity anomaly extraction module, and fusion reconstruction module;

[0272] The gravity measurement module is used to obtain ground gravity measurement data and calculate the original ground gravity anomaly of its discrete points on the ground;

[0273] The long-wave gravity anomaly extraction module is used to introduce the global gravity field model and calculate the long-wave gravity anomaly components of discrete points on the ground in real time;

[0274] The shortwave gravity anomaly extraction module is used to obtain a terrain data sample set; a deep fully connected neural network (FC-DNN) is introduced; the sample locations and terrain features in the terrain data sample set are used as input, and the shortwave gravity anomaly therein is used as a supervision target to train the FC-DNN to learn the mapping relationship between the terrain features and the shortwave gravity anomaly; the trained FC-DNN is used as a shortwave gravity field fast solution model, and the shortwave gravity anomaly components are obtained based on the terrain features of discrete points on the ground in real time.

[0275] The residual gravity anomaly extraction module is used to remove the long-wave and short-wave gravity anomaly components of the ground discrete points from their original ground gravity anomalies to obtain their residual gravity anomalies; construct a multi-source data sample set, train a CNN to model the nonlinear relationship between the multi-source data and the residual gravity anomaly; use the trained CNN as a rapid solution model for the residual gravity anomaly, and obtain its residual gravity anomaly components based on real-time prediction of the multi-source data of the ground discrete points;

[0276] The fusion reconstruction module is used to fuse the residual gravity anomaly component, the short-wave gravity anomaly component, and the long-wave gravity anomaly component of the ground discrete points to complete the modeling of the regional ground gravity field model.

[0277] Obviously, the above embodiments of the present invention are merely examples for the purpose of illustrating the present invention, and are not intended to limit the embodiments of the present invention. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention shall be included within the scope of protection of the claims of the present invention.

Claims

1. A real-time regional gravity field modeling method based on deep neural network, characterized by: The specific steps include: Obtain ground gravity measurement data and calculate the original ground gravity anomaly of its discrete points; The global gravity field model is introduced to calculate the long-wave gravity anomaly components of discrete points on the ground in real time; A terrain data sample set is obtained; a deep fully connected neural network (FC-DNN) is introduced; the sample locations and terrain features in the terrain data sample set are used as input, and the shortwave gravity anomaly is used as the supervision target to train the FC-DNN to learn the mapping relationship between terrain features and shortwave gravity anomaly. The trained FC-DNN is used as a fast shortwave gravity field solution model, and the shortwave gravity anomaly component is obtained in real time based on the terrain features of discrete ground points. The long-wave and short-wave gravity anomaly components of the ground discrete points are removed from the original ground gravity anomaly to obtain the residual gravity anomaly. Constructing a multi-source data sample set and training a CNN to model the nonlinear relationship between the multi-source data and the residual gravity anomaly; The trained CNN is used as a fast solution model for residual gravity anomaly, and its residual gravity anomaly components are obtained by real-time prediction based on multi-source data of discrete ground points. The residual gravity anomaly component of discrete ground points, the short-wave gravity anomaly component, and the long-wave gravity anomaly component are integrated to complete the modeling of the regional ground gravity field model.

2. The real-time regional gravity field modeling method based on deep neural network according to claim 1 is characterized in that: The original ground gravity anomaly is calculated. The specific steps are: Obtain ground gravity measurement data for each sample, including longitude, latitude, surface elevation, height of discrete ground points, and the difference between the height of discrete ground points and the surface elevation; Perform preprocessing to eliminate systematic errors and outliers in ground gravity measurement data; The normal gravity is subtracted from the preprocessed gravity measurement data to obtain the original ground gravity anomaly of its discrete points.

3. The real-time regional gravity field modeling method based on deep neural network according to claim 1 is characterized in that: The global gravity field model is introduced to construct the long-wave gravity anomaly components of discrete points on the ground. The specific steps are as follows: Selecting the global gravity field model, the gravitational potential V at a point P outside the Earth is defined by its geocentric distance r, geocentric colatitude θ and longitude λ: Where GM is the Earth's gravitational constant, a is the fully normalized spherical harmonic coefficient Related scale factors, in addition: in, is a fully normalized Legendre function of the first kind, with order n and degree |m|. The first-order term (n=1) in the spherical harmonic expansion represents the offset of the center of gravity. The disturbance potential T at point P(r,θ,λ) is defined as the difference between the actual gravity potential of the Earth and the normal gravity potential; The perturbation potential is expressed by spherical harmonic expansion: The zero-order term in the formula reflects the gravitational effect produced by the total mass of the Earth. If the mass of the reference ellipsoid is equal to the actual mass of the Earth, this term is zero. The surface free-air gravity anomaly Δg is defined as the actual gravity acceleration at point P minus the normal gravity acceleration at the corresponding geoid point Q: Δg=|g P |-|c Q | Assuming that the gravitational potential on the geoid is equal to the normal gravitational potential on the reference ellipsoid, the basic relationship between gravity anomaly and disturbance potential is obtained: in represents the radial derivative of the perturbation potential, is a correction term related to the distance, and the three ε terms represent three types of error corrections; Let Δg c Represents the gravity anomaly after all systematic corrections have been applied: get: The long-wave gravity anomaly component is further obtained: in From this we get The spherical harmonic expansion expression of , whose infinite series in the external space of the Brillouin sphere has convergence.

4. The real-time regional gravity field modeling method based on deep neural network according to claim 1 is characterized in that: Obtain a sample set of terrain data. The specific steps are as follows: Define source quality models; collect high-resolution terrain feature data and extract terrain features including surface elevation data; The high-resolution terrain feature data is filtered to obtain reference terrain data, and the difference between the reference terrain data and the high-resolution terrain feature data is used as the residual terrain model (RTM); Randomly select discrete points from the terrain data, record the sample position of each discrete point including longitude, latitude and altitude, and calculate its shortwave gravity anomaly; The RTM mass distribution is divided into several regions according to the distance from the calculation point, and the shortwave gravity anomaly is approximately calculated for each region using geometric bodies; The terrain data sample set is obtained by fusing terrain features, shortwave gravity anomalies, and sample locations.

5. The real-time regional gravity field modeling method based on deep neural network according to claim 1 is characterized in that: The FC-DNN includes an input layer, L hidden layers, and an output layer. During training, the FC-DNN learns the mapping relationship between terrain features and shortwave gravity anomalies in an end-to-end manner, specifically: Assume that there are N samples of terrain data and adopt batch training method. Set the number of samples in each batch to k. According to N / k=q, divide the data into q batches. Each batch performs forward operation to obtain the loss function, and then backpropagates to update the weights. This process is called one iteration. The forward propagation of each layer is defined as: H (0) =X Z (l) =W (l) H (l-1) +B (l) H (l) =σ(Z (l) ) Among them, X is the matrix of input samples, with dimension [m×k], W (l) is the weight, B (l) is the bias, which is broadcast to k columns during calculation to match the batch size; The batch loss function for training is: Where Y is the true label matrix with dimension [n L ×k],H (L) Is the output matrix of the output layer, with dimension [n L ×k],||·|| F represents the Frobenius norm; During the backpropagation process, the error term of the output layer is calculated as: Where ⊙ represents the element-by-element multiplication of the matrix; For the hidden layer (l = L-1, L-2, ..., 1), the recursive formula of the error term is: Δ (l) =(W (l+1) ) T Δ (l+2) ⊙σ′(Z (l) ) The batch gradient is calculated as: Among them 1 k is a k-dimensional all-one column vector used to sum the gradients over the batch dimension; The parameter update formula is: In batch training, all batches undergo one iteration and complete one cycle, and training continues until the set threshold is reached to obtain the trained FC-DNN.

6. The real-time regional gravity field modeling method based on deep neural network according to claim 5 is characterized in that: The FC-DNN training also introduces a terrain regularization method to further improve the generalization ability and prediction accuracy of the network; the loss function of the regularization method is: Where a is a trainable array with a dimension of [1×k]. The mean and variance used to initialize it are the mean and variance of the gravity anomaly and terrain ratio in the training area. Predict gravity anomalies for the training area, h is H (L) The corresponding surface elevation value.

7. The real-time regional gravity field modeling method based on deep neural network according to claim 1 is characterized in that: The structure of CNN includes: Feature extraction part: includes a convolution layer, batch normalization layer, activation function layer, and pooling layer; Position information fusion part: flatten the feature map output by the feature extraction part and then splice it with the position information of the discrete points on the ground; Regression prediction part: includes at least one fully connected layer, which is used to map the fused feature vector into the residual gravity anomaly prediction value.

8. The real-time regional gravity field modeling method based on deep neural network according to claim 7 is characterized in that: The specific training process of CNN is as follows: Assume that the input of the lth convolutional layer in the feature extraction part is h (l-1) , the convolution kernel is K (l) , bias is b (l) , define the output of the convolutional layer as: h (l) =f(Conv(h (l-1) ,K (l) )+b(l)) Where f(·) is the activation function and Conv(·) is the convolution operation, which is defined as: Among them O i,j is the value of the output feature map at position (i, j), k is the convolution kernel size, and b is the bias term; The batch normalization layer is defined as: Among them, μ B and is the mean and variance of the batch, γ and β are learnable scaling and translation parameters, is a constant; The output of the batch normalization layer passes through the Tanh activation function layer and then the pooling layer to further reduce the feature map size by downsampling: Where R is the pooling window, and features are extracted layer by layer by stacking multiple convolutional layers and pooling layers; In the position information fusion part, the extracted feature map is flattened and then spliced ​​with the position of the discrete points on the ground to form a complete feature vector; The complete feature vector is flattened into a one-dimensional vector, and the fully connected layer is connected to the output layer for regression task; Introducing the weighted loss function: Construct a gradient descent method to minimize the loss function E: Among them, w i is the error weight, according to y i The distribution definition of y i When the distribution is uneven, values ​​in intervals with fewer numerical distributions are assigned greater weights; Reshape the error into the shape of the feature map: For the max pooling layer, the gradient is only passed to the location of the maximum value in the pooling window: Where |R| is the size of the pooling window; The gradient of the batch normalization layer is calculated as: The gradient of the loss function to the convolution kernel in the output gradient of the convolution layer l is: The gradient with respect to the bias is: The gradient of the input feature map is: Calculated As the error input of the previous layer, continue backpropagation training; Update the network parameters using gradient descent: The input multi-source data sample set is divided into a training set and a test set. The training set is divided into multiple batches and trained for several cycles. Some samples are selected from the training set as a validation set. The network accuracy is verified after each cycle of training. At the same time, the hyperparameters are adjusted to ensure that the accuracy meets expectations and there is no overfitting. The trained CNN is obtained.

9. The real-time regional gravity field modeling method based on deep neural network according to claim 8, characterized in that: In order to solve the problem of gradient disappearance and network degradation faced by CNN, at least one residual block is introduced in the feature extraction part, and its structure is as follows: y=σ(F(x,{W i })+x) Where σ represents the nonlinear transformation composed of the convolution layer, batch normalization layer and activation function; x is the feature map extracted by the feature extraction part, and y is the complete feature map output by the residual block; For the residual block, the error backpropagation is expressed as: h (l) =f(Conv(h (l-1) ,K (l) )+b(l))。 10. A real-time regional gravity field modeling system based on deep neural networks, characterized by: It includes gravity measurement module, long-wave gravity anomaly extraction module, short-wave gravity anomaly extraction module, residual gravity anomaly extraction module, and fusion reconstruction module; The gravity measurement module is used to obtain ground gravity measurement data and calculate the original ground gravity anomaly of its discrete points on the ground; The long-wave gravity anomaly extraction module is used to introduce the global gravity field model and calculate the long-wave gravity anomaly components of discrete points on the ground in real time; The shortwave gravity anomaly extraction module is used to obtain a terrain data sample set; a deep fully connected neural network (FC-DNN) is introduced; the sample locations and terrain features in the terrain data sample set are used as input, and the shortwave gravity anomaly therein is used as a supervision target to train the FC-DNN to learn the mapping relationship between the terrain features and the shortwave gravity anomaly; the trained FC-DNN is used as a shortwave gravity field fast solution model, and the shortwave gravity anomaly components are obtained based on the terrain features of discrete points on the ground in real time. The residual gravity anomaly extraction module is used to remove the long-wave and short-wave gravity anomaly components of the ground discrete points from the original ground gravity anomaly to obtain the residual gravity anomaly; Constructing a multi-source data sample set and training a CNN to model the nonlinear relationship between the multi-source data and the residual gravity anomaly; The trained CNN is used as a fast solution model for residual gravity anomaly, and its residual gravity anomaly components are obtained by real-time prediction based on multi-source data of discrete ground points. The fusion reconstruction module is used to fuse the residual gravity anomaly component, the short-wave gravity anomaly component, and the long-wave gravity anomaly component of the ground discrete points to complete the modeling of the regional ground gravity field model.

Images