A time-domain airborne electromagnetic three-dimensional inversion method based on POD reduction

Abstract

The application belongs to the technical field of geophysical inversion, and particularly relates to a time domain airborne electromagnetic three-dimensional inversion method based on POD reduction, which comprises the following steps: calculating three-dimensional electromagnetic responses corresponding to each undersampling transmitting source, iteratively reconstructing complete time domain prediction responses corresponding to each transmitting source position, performing full-order Gauss-Newton inversion to obtain a POD source resistivity model, calculating POD bases and mapping to an internal inversion parameter space to obtain a reduction base, using the reduction base to obtain a reduction model update direction, generating a reduction trial resistivity model through reduction line search, and iteratively accepting the reduction trial resistivity model meeting a target function value drop condition as a current inversion model. The method can improve the numerical stability of the inversion process while improving the calculation efficiency of the three-dimensional time domain airborne electromagnetic inversion, and can maintain the reliability of the inversion result.

Description

Technical Field

[0001] This application belongs to the field of geophysical inversion technology, specifically a time-domain airborne electromagnetic three-dimensional inversion method based on POD order reduction. Background Technology

[0002] The time-domain airborne electromagnetic method uses a flight platform equipped with transmitting and receiving devices to measure the attenuation response of the secondary field generated by the underground medium after the transmitting current is turned off. It features high coverage efficiency and sensitivity to low-resistivity anomalies. Three-dimensional inversion can fully consider complex geological bodies and three-dimensional electromagnetic coupling effects, but each iteration usually requires multiple large-scale forward models, sensitivity calculations, and solutions to high-dimensional linear equations, resulting in large computational and storage requirements.

[0003] When the survey area contains a large number of emission source locations, performing 3D forward modeling on all emission sources in every model evaluation would result in repeated line search and inversion iterations, leading to high computational costs. Airborne electromagnetic responses typically exhibit spatial correlation along the survey line or area. Therefore, while maintaining a representative spatial distribution, a subset of emission sources can be selected for calculation, and the compressibility of the response in the sparse transform domain can be used to recover the predicted response at uncalculated locations. However, if compressed sampling, response reconstruction, full-order inversion, and model order reduction are independent of each other, inconsistencies in data arrangement, sampling locations, and objective function evaluation can easily arise.

[0004] Gauss-Newton inversions often use data covariance to weight the residuals and constrain model updates using model covariance or smoothing operators. To improve numerical conditions, resistivity model parameters can be transformed to an internal inversion parameter space. In this case, the updates in the forward model space and the updates in the internal parameter space are coupled through the square root of the model covariance operator. If the POD basis obtained on a regular resistivity grid is directly used for internal parameter updates, inconsistencies in grid dimension, cell arrangement, and parameter meanings will arise simultaneously.

[0005] On the other hand, time-domain airborne electromagnetic sensitivity is often stored as a local window based on the transmitter source, while the model parameters involved in the inversion are stored in a global internal cell array. If the local sensitivity is directly flattened into a global vector, or if the row order of transmitter source, time step, measurement point, component, and time channel is not unified, it is easy to cause out-of-bounds errors, misaligned projections, and abnormal changes in the objective function. Summary of the Invention

[0006] This application provides a time-domain airborne electromagnetic three-dimensional inversion method based on POD order reduction, which solves the problems of inconsistency between the compressed sampling response and the complete observation data, inconsistency between the POD regular grid model and the internal inversion parameter dimension, loss of model covariance mapping, and inability of order reduction update to maintain the original objective function constraints.

[0007] This application provides a time-domain airborne electromagnetic three-dimensional inversion method based on POD-reduced order, including: Determine the set of undersampled emission sources based on the spatial location of the emission sources in the test area; A local computational region is established for each undersampled emission source in the set of undersampled emission sources, and the three-dimensional electromagnetic response corresponding to each undersampled emission source is calculated. The three-dimensional electromagnetic response is converted into an undersampled time-domain predicted response for each measurement time channel, and the complete time-domain predicted response corresponding to each emission source position is iteratively reconstructed in the sparse transform domain. A full-order Gaussian-Newton inversion is performed. In each iteration, the full-order weighted residual and full-order weighted sensitivity matrix are calculated based on the predicted response in the complete time domain. The gradient of the objective function and the Gaussian-Newton approximation Hessian matrix in the full-order Gaussian-Newton inversion are constructed based on the full-order weighted residual and full-order weighted sensitivity matrix, and the full-order Gaussian-Newton equation is formed. The resistivity model update direction is obtained by solving the full-order Gaussian-Newton equation. The accepted full-order trial resistivity model is determined by full-order line search. When the preset order reduction switching condition is met, the last accepted full-order trial resistivity model before order reduction is selected as the POD source resistivity model. Multiple snapshots are constructed centered on the POD source resistivity model, forming a snapshot matrix. Singular value decomposition is then performed on the snapshot matrix to obtain the POD basis. The POD basis is mapped to the internal inversion parameter space to obtain a reduced-order basis; Using the reduced-order basis, the full-order weighted residual and the full-order weighted sensitivity matrix are projected onto the reduced-order space, and the reduced-order Gauss-Newton equation is solved to obtain the update direction of the reduced-order model; The order reduction model update direction is generated by searching the order reduction line to generate an order reduction trial resistivity model, and the complete time domain predicted response of the order reduction trial resistivity model is calculated; only the order reduction trial resistivity model that satisfies the objective function value descent condition is accepted as the current inversion model; When the updated current inversion model meets the preset termination condition, the final inversion model is output.

[0008] Furthermore, based on the spatial location of the emission sources in the survey area, a set of undersampled emission sources is determined, including: The Poisson disk sampling method is used to select the emission source positions in the two-dimensional spatial grid of the test area, such that the spatial distance between any two selected emission source positions is not less than a preset minimum spacing threshold, thus obtaining the undersampled emission source set.

[0009] Furthermore, a local computational region is established for each undersampled emission source in the set of undersampled emission sources, and the three-dimensional electromagnetic response corresponding to each undersampled emission source is calculated, including: Centered on the current undersampled emission source location, the local effective influence area is determined based on the skin depth; A local grid is established within the local effective influence area, and an extended cell layer is provided outward from the boundary of the local grid to absorb boundary reflections; Solving the secondary electric field diffusion equation yields the secondary electric field distribution at each node on the local grid. This distribution is then superimposed with the background field to obtain the three-dimensional electromagnetic response of the current undersampled emission source.

[0010] Further, the three-dimensional electromagnetic response is converted into an undersampled time-domain predicted response for each measurement time channel, and the complete time-domain predicted response corresponding to each emission source position is iteratively reconstructed in the sparse transform domain, including: The frequency domain response of each undersampled emission source is transformed into the frequency-time domain to obtain the undersampled time domain predicted response of each measurement time channel; For each measurement time channel, the undersampled time domain predicted response is arranged into a two-dimensional response matrix according to the spatial location of the emission source, the unsampled positions are set to zero, and a binary sampling position indication matrix of the same size as the two-dimensional response matrix is ​​established. Using curvelet transform as a sparse transform, the two-dimensional response matrix of each measurement time channel is iteratively reconstructed. Based on the binary sampling position indication matrix, the response values ​​of the sampled positions in the two-dimensional response matrix are forcibly replaced with the corresponding values ​​in the original undersampled time domain predicted response, thus completing the data reconstruction of each measurement time channel. After iterative convergence, the data reconstruction results are assembled in a predetermined order into complete time-domain predicted responses corresponding to each emission source location.

[0011] Furthermore, in each iteration, the full-order weighted residual and full-order weighted sensitivity matrix are calculated based on the predicted response in the complete time domain. The gradient of the objective function and the Gauss-Newton approximation Hessian matrix in the full-order Gauss-Newton inversion are constructed based on these matrixes, forming the full-order Gauss-Newton equation. Solving the full-order Gauss-Newton equation yields the resistivity model update direction, including: The full-order weighted residuals are calculated based on the difference between the observed data and the predicted response in the complete time domain, combined with the data covariance matrix. Calculate the full-order weighted sensitivity of the predicted response in the full time domain to the resistivity model parameters under the current resistivity model; when the full-order weighted sensitivity is stored in the local window of the emission source, extract the corresponding row from the global network model space through index mapping, and combine them into a full-order weighted sensitivity matrix in the order of emission source, time step, measurement point, component and time channel; The gradient of the objective function is calculated from the full-order weighted residual and the full-order weighted sensitivity matrix. The Gauss-Newton approximation Hessian matrix is ​​constructed by multiplying the transpose of the full-order weighted sensitivity matrix with the full-order weighted sensitivity matrix, plus the product of the regularization parameter and the data covariance matrix. The linear equation system of the Gauss-Newton approximation Hessian matrix and the gradient of the objective function is solved to obtain the resistivity model update direction.

[0012] Furthermore, multiple snapshots are constructed centered on the POD source resistivity model to form a snapshot matrix. Singular value decomposition is then performed on the snapshot matrix to obtain the POD basis, including: The spatial regions in the POD source resistivity model where the resistivity value deviates from the background resistivity by a preset value are marked as abnormal regions. Using the abnormal regions as perturbation objects, a set of snapshots is generated using multiple perturbation modes. The generated snapshots are arranged in columns to form a snapshot matrix. The snapshot matrix is ​​then processed to remove the mean and center the snapshots so that the mean of each snapshot is zero. Singular value decomposition is performed on the mean-removed snapshot matrix. The resulting singular values ​​are arranged in descending order. The square of each singular value represents the energy share of the corresponding mode in the snapshot matrix. The energy share of the first few singular values ​​is accumulated sequentially. When the proportion of the accumulated energy to the total energy first reaches the preset holding threshold, the number of modes contained in the current accumulation is used as the reduced dimension. Extract the left singular vectors corresponding to the selected reduced dimension from the singular value decomposition results to form the POD basis.

[0013] Further, the POD basis is mapped to the internal inversion parameter space to obtain a reduced-order basis, including: Based on the coordinate correspondence between regular grid cells and internal inversion cells, a constraint mapping from the regular grid index to the internal inversion cell index is established. The constraint mapping is used to filter and rearrange the POD base, remove redundant rows corresponding to non-inversion cells, and rearrange the rows corresponding to internal inversion cells according to the storage order of the internal inversion cells to obtain the POD base in the internal cell space. Using the square root of model covariance as the mapping operator, the mapping equation is solved for each column of the POD basis in the internal cell space to obtain the reduced-order basis of the internal inversion parameter space.

[0014] Further, the basis mapping consistency error is calculated. The basis mapping consistency error is the relative deviation between the POD basis in the internal cell space and the reduced-order basis result of the model covariance square root operator acting on the internal inversion parameter space. When the basis mapping consistency error exceeds a preset threshold, the reduced-order iteration is stopped and the iteration is rolled back to the full-order iteration.

[0015] Furthermore, using the reduced-order basis, the full-order weighted residuals and full-order weighted sensitivity matrices are projected onto the reduced-order space, and the reduced-order Gauss-Newton equations are solved to obtain the updated direction of the reduced-order model, including: Multiply the full-order weighted sensitivity matrix by the reduced-order basis to obtain the reduced-order weighted sensitivity matrix; Multiplying the reduced-order weighted sensitivity matrix by the full-order weighted residual yields the reduced-order negative gradient vector; The left-hand coefficient matrix of the reduced-order Gauss-Newton equation is formed by multiplying the transpose of the reduced-order weighted sensitivity matrix with the product of the regularization parameter and the transpose of the reduced-order basis with the product of the reduced-order basis; the reduced-order Gauss-Newton equation is constructed by using the reduced-order negative gradient vector as the right-hand term of the reduced-order equation. Solve the reduced-order Gauss-Newton equation to obtain the update direction of the reduced-order model.

[0016] Furthermore, a reduced-order trial resistivity model is generated by searching for the reduced-order line in the updated direction of the reduced-order model, and the complete time-domain predicted response of the reduced-order trial resistivity model is calculated; only reduced-order trial resistivity models that satisfy the objective function value descent condition are accepted as the current inversion model, including: With a step size of 1 as the initial value, a trial reduction coefficient is generated along the update direction of the reduction model. Based on the linear combination of the current internal inversion model, the reduction basis, and the trial reduction coefficient, a trial internal inversion model is generated. Based on the mapping relationship between the internal inversion parameters and the resistivity model parameters, the trial internal inversion model is converted into a reduced-order trial resistivity model. Calculate the complete time-domain predicted response of the reduced-order trial resistivity model; The objective function value is calculated based on the difference between the observed data and the predicted response in the complete time domain, and it is determined whether the objective function value satisfies the descent condition relative to the objective function value of the current inversion model. If the descent condition is met, the trial resistivity model is accepted as the current inversion model, and the corresponding full time domain prediction response, objective function value, and normalized root mean square error are updated. If the descent condition is not met, the trial order reduction coefficient, trial internal inversion model, and reduced-order trial resistivity model are regenerated after reducing the step size by a preset ratio. The full time domain prediction response is recalculated and re-evaluated until the descent condition is met or the maximum number of attempts is reached. If no reduced resistivity model satisfies the descent condition after reaching the maximum number of attempts, the current inversion model remains unchanged, and the search for the reduced line is terminated.

[0017] Compared with the prior art, this application has at least the following beneficial effects: This application reduces the number of 3D forward modeling calculations required for each model evaluation by spatially undersampling the emission source location and reconstructing the time-domain predicted response of the complete survey area in the sparse transform domain. It constructs a snapshot matrix and extracts the POD basis based on the last accepted full-order trial resistivity model before order reduction, enabling the reduced-order space to reflect the main structural features of the current inversion model, reducing the degrees of freedom for updating the Gauss-Newton model, and decreasing the computational cost of solving the linear equations. Furthermore, it establishes a consistent mapping between the internal cell space of the regular grid, the internal inversion cell space, and the internal inversion parameter space. This approach avoids mismatches in dimensionality, index order, and physical meaning between different model spaces, while maintaining consistency between the reduced-order model update and the model covariance constraint in the original objective function. Furthermore, the same emission source sampling index and response reconstruction process are used in full-order iteration, POD reduced-order iteration, and line search trial resistivity model evaluation. Verification is performed by combining response reconstruction error, basis mapping relative deviation, reduced-order linear solution residuals, and the original objective function descent condition. This improves the computational efficiency of three-dimensional time-domain airborne electromagnetic inversion while enhancing the numerical stability of the inversion process and maintaining the reliability of the inversion results. Attached Figure Description

[0018] Figure 1 A flowchart illustrating the overall process of the time-domain airborne electromagnetic three-dimensional inversion method based on POD order reduction provided in this application embodiment; Figure 2 This is a schematic diagram of the co-mapping between the POD basis of the regular grid and the internal inversion parameter space provided in the embodiments of this application, wherein (a) is the POD basis of the regular grid, (b) is the result of effective cell selection and parameter order rearrangement of the POD basis, (c) is the POD basis under the rearranged internal cell space, (d) is the reduced-order basis of the internal inversion parameter space, (e) is the full-order weighted sensitivity matrix, and (f) is the reduced-order weighted sensitivity matrix. Figure 3 This is a schematic diagram of a real resistivity model provided for an embodiment of this application; Figure 4 A schematic diagram of the inversion model obtained by the method of this application as provided in an embodiment of this application. Detailed Implementation

[0019] To make the objectives, technical solutions, and advantages of this application clearer, the following detailed description is provided in conjunction with embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the scope of this application.

[0020] In this embodiment of the application, the data preparation stage of the entire inversion method includes: acquiring time-domain airborne electromagnetic observation data, a priori resistivity model and three-dimensional mesh, setting data error, model covariance, regularization parameter, full-order iteration number, spatial sampling parameter and inversion termination condition, and establishing the model transformation relationship between resistivity model parameters and Gauss-Newton inversion parameters.

[0021] The observation data consists of secondary field attenuation response data obtained by the time-domain airborne electromagnetic system during flight measurements in the actual survey area. Each observation data corresponds to a transmitter location, a receiver location, a time channel, and a field component. The observation data are the known conditions and fitting target for the inversion problem. The essence of the inversion is to infer the distribution of underground resistivity from the electromagnetic response measured on the surface. Without observation data, the inversion loses its solution objective, and the entire iterative process cannot proceed.

[0022] The prior resistivity model is the initial inversion model, typically a uniform half-space model or a layered medium model determined by other prior geological information. In this embodiment, the prior model is a 100Ω·m uniform half-space model. The three-dimensional mesh is a regular hexahedral mesh covering the entire survey area, used for discretization of the electromagnetic field in forward modeling calculations.

[0023] A three-dimensional mesh serves as a numerical computational medium for discretizing underground space. The electromagnetic field diffusion equations are solved discretized on mesh cells. Without a three-dimensional mesh, it is impossible to define the spatial location and dimensions of model parameters, and numerical forward modeling is also impossible. The three-dimensional mesh is a regular hexahedral mesh, used for discretizing the electromagnetic field in forward modeling. It also defines a complete model parameter space. Subsequent selection of internal inversion cells, snapshot construction, and POD basis extraction are all completed on this three-dimensional mesh.

[0024] The data covariance matrix is ​​used to weight the observations, with its diagonal elements representing estimates of the variance of each observation. Earlier, higher-quality data are assigned smaller standard deviations (higher weights), while later, lower-quality data are assigned larger standard deviations (lower weights), ensuring the inversion focuses more on reliable data. Through covariance weighting, data with different dimensions and noise levels are made comparable within the objective function.

[0025] The model covariance matrix describes the spatial correlation between model parameters and the smoothness of parameter variations. Specifically, it is a Gaussian covariance function. The correlation length and model variance together determine the spatial variation characteristics of resistivity. The square root operator of the model covariance matrix is ​​used to establish the transformation relationship between resistivity model parameters and internal inversion parameters.

[0026] Regularization parameters are used to balance the weights between data fitting terms and model constraint terms in the objective function. The larger the regularization parameter, the stronger the model constraint and the smoother the inversion result; the smaller the regularization parameter, the more dominant the data fitting becomes, and the model may produce unnecessary jitter.

[0027] Set the number of full-order iterations. The number of full-order iterations refers to the number of full-order Gaussian-Newton iterations performed before entering the POD reduction inversion. Its value should ensure that the model converges sufficiently to the vicinity of the true solution, so that the POD basis can be constructed based on a relatively reliable source model, while avoiding excessive computational overhead.

[0028] Spatial sampling parameters are set. These parameters control the sampling density and minimum spacing constraint of the undersampled emitter set. The spatial sampling parameters directly affect the balance between forward modeling computation and response reconstruction accuracy. A higher sampling ratio results in greater computational complexity but higher reconstruction accuracy; a lower sampling ratio leads to more significant improvements in computational efficiency but may increase reconstruction error. In this embodiment, the sampling spacing control parameter is set to 1.09, corresponding to a sampling ratio of approximately 41%.

[0029] Set inversion termination conditions. These conditions control when the iteration stops, preventing infinite loops or wasted computational resources. They include at least one of the following: the normalized root mean square error reaches a target value; the relative change in the objective function is below a preset threshold; the reduced negative gradient norm is below a preset threshold; the model update amount is below a preset threshold; or the maximum number of iterations is reached. When any termination condition is met, the inversion iteration stops and the current inversion model is output.

[0030] The model transformation relationship is defined in the internal inversion parameter space as follows: ; ; ; in, For three-dimensional time-domain airborne electromagnetic forward modeling operators, For the data covariance matrix, The model covariance matrix, For regularization parameters, The objective function value, For data fitting terms, For model constraint terms, For observation data, These are resistivity model parameters. This is the a priori resistivity model, where T is the transpose.

[0031] After performing model covariance transformation, internal inversion parameters are introduced. : , , Among them, internal inversion parameters With resistivity model parameters satisfy: ,in For the prior circuit rate model, The square root operator for the model covariance; The inverse operator of the square root of the model covariance operator; Substituting the above equation into the objective function, we obtain the objective function in the internal inversion parameter space: , The objective function includes a data fit term weighted by the data covariance. and model constraint terms controlled by regularization parameters. ; Let be the inverse matrix of the square root of the data covariance.

[0032] Based on the data prepared above, this application provides a time-domain airborne electromagnetic three-dimensional inversion method based on POD-reduced order, including: S101, Determine the set of undersampled emission sources based on the spatial location of the emission sources in the test area; S102, establish a local computational region for each undersampled emission source in the set of undersampled emission sources, and calculate the three-dimensional electromagnetic response corresponding to each undersampled emission source; S103, the three-dimensional electromagnetic response is converted into an undersampled time-domain prediction response for each measurement time channel, and the complete time-domain prediction response corresponding to each emission source position is iteratively reconstructed in the sparse transform domain. S104, perform full-order Gauss-Newton inversion. In each iteration, calculate the full-order weighted residual and full-order weighted sensitivity matrix based on the full-time domain predicted response. Construct the gradient of the objective function and the Gauss-Newton approximation Hessian matrix in the full-order Gauss-Newton inversion based on the full-order weighted residual and full-order weighted sensitivity matrix, and form the full-order Gauss-Newton equation. Obtain the resistivity model update direction by solving the full-order Gauss-Newton equation. Determine the accepted full-order trial resistivity model through full-order line search. When the preset order reduction switching condition is met, select the last accepted full-order trial resistivity model before order reduction as the POD (Intrinsic Orthogonal Decomposition) source resistivity model. S105, construct multiple snapshots centered on the POD source resistivity model and form a snapshot matrix. Perform singular value decomposition on the snapshot matrix to obtain the POD basis. S106, Map the POD basis to the internal inversion parameter space to obtain a reduced-order basis; S107. Using the reduced-order basis, project the full-order weighted residual and the full-order weighted sensitivity matrix onto the reduced-order space, and solve the reduced-order Gauss-Newton equation to obtain the updated direction of the reduced-order model. S108, the order reduction model update direction is generated by order reduction line search to generate an order reduction trial resistivity model, and the complete time domain predicted response of the order reduction trial resistivity model is calculated; only the order reduction trial resistivity model that satisfies the objective function value descent condition is accepted as the current inversion model; S109. When the updated current inversion model meets the preset termination condition, output the final inversion model. In step S101 of this application, determining the set of undersampled emission sources based on the spatial location of the emission sources in the test area includes: selecting emission source locations in a two-dimensional spatial grid of the test area using the Poisson disk sampling method, such that the spatial distance between any two selected emission source locations is not less than a preset minimum spacing threshold, thereby obtaining the set of undersampled emission sources.

[0033] During flight in the survey area, the time-domain airborne electromagnetic observation system collects data along the survey line at close transmission intervals. Taking one embodiment as an example, the observation system includes 861 transmitter locations, evenly distributed along the survey line direction (x-direction) and perpendicular to the survey line direction (y-direction). If three-dimensional forward modeling is performed on all 861 transmitters in each model evaluation, the computational load is extremely high, resulting in low inversion iteration efficiency.

[0034] However, the airborne electromagnetic response exhibits significant correlation across the spatial range of the survey area. The electromagnetic response changes from adjacent emission sources are continuous and gradual, without drastic abrupt changes. This indicates that the electromagnetic response data possesses redundancy in the spatial domain, and the complete response of the entire survey area can be reconstructed from the responses of a small number of sampling points. Based on compressed sensing theory, if the signal exhibits sparsity in a certain transform domain, the complete signal can be recovered with high precision from observations at a sampling rate far lower than the Nyquist sampling rate.

[0035] Therefore, this application adopts a spatial undersampling strategy: only three-dimensional forward modeling is performed on some of the emission source locations in the test area, and then the complete time domain prediction response of all emission sources is obtained through subsequent sparse transform domain reconstruction.

[0036] In this embodiment, the Poisson disk sampling method is used to determine the set of undersampled emission sources. Poisson disk sampling is a random sampling method whose core feature is that the spatial distance between any two sampling points is not less than a preset minimum spacing threshold, ensuring that the sampling points are evenly distributed in space and do not cluster together.

[0037] The specific sampling process is as follows: The two-dimensional spatial plane of the survey area is divided into a regular square grid to form a two-dimensional spatial network. The grid side length is set as the minimum spacing threshold divided by 1 / 2. Using this grid size ensures that if two sampling points fall into adjacent grids, the distance between them will always be less than the minimum spacing threshold. Therefore, only checking adjacent grids is needed to determine if a candidate point meets the spacing constraint, significantly improving sampling efficiency. In this embodiment, the test area has 861 transmitters at 21 positions in the x-direction and 41 positions in the y-direction, forming a two-dimensional grid.

[0038] A first sampling point is randomly generated within the two-dimensional space of the survey area as a seed point, and it is added to both the processing array and the result array. The processing array stores the sampling points that need to be expanded, and the result array stores the final set of undersampled emission sources.

[0039] When the processing array is not empty, a sampling point is randomly selected and removed from the processing array. Centered on this point, a predetermined number of candidate sampling points are randomly generated within a ring-shaped region with an inner diameter equal to the minimum spacing threshold and an outer diameter equal to twice the minimum spacing threshold. Each candidate point undergoes a validity check: it is checked whether the point is within the spatial range of the survey area; and it is checked whether the distance between the point and all existing sampling points in the result array is not less than the minimum spacing threshold. When checking the distance to existing sampling points, the grid structure only needs to check points in the surrounding adjacent grids, eliminating the need to traverse all existing points, thus significantly improving computational efficiency.

[0040] If a candidate point passes the validity check, it is added to both the processing array and the result array. If no valid candidate point can be generated after a preset number of random attempts, the expansion of that point is abandoned. The above process is repeated until the processing array is empty. At this point, the set of points in the result array is the set of undersampled emission sources.

[0041] After selecting the spatial locations of the undersampled emission source set, it is necessary to establish an index mapping relationship between this undersampled emission source set and the complete emission source sequence. Specifically: The transmitters in the survey area are arranged in a two-dimensional index matrix according to their spatial location. Each transmitter has a two-dimensional coordinate system, including a row number and a column number, where the row number corresponds to a position in one horizontal direction and the column number corresponds to a position in another horizontal direction. This two-dimensional coordinate system is then converted into a one-dimensional sequence number in the complete transmitter sequence.

[0042] Construct a logical array of length equal to the total number of complete emitters. The array elements are Boolean: if the corresponding emitter is selected into the undersampled emitter set, mark it as true; otherwise, mark it as false. The resulting undersampled emitter set remains consistent in the full-order stage, POD stage, and line search model evaluation of the same inversion.

[0043] In step S102 of this application, a local computational region is established for each undersampled emission source in the set of undersampled emission sources, and the three-dimensional electromagnetic response corresponding to each undersampled emission source is calculated, including: Centered on the current undersampled emission source location, the local effective influence area is determined based on the skin depth; A local grid is established within the local effective influence area, and an extended cell layer is provided outward from the boundary of the local grid to absorb boundary reflections; Solving the secondary electric field diffusion equation yields the secondary electric field distribution at each node on the local grid. This distribution is then superimposed with the background field to obtain the three-dimensional electromagnetic response of the current undersampled emission source.

[0044] Time-domain airborne electromagnetic three-dimensional forward modeling is the process of obtaining the electromagnetic response at each receiving point by numerically solving Maxwell's equations under given resistivity model and transmitter parameters. Traditional three-dimensional forward modeling is usually solved on a one-time grid covering the entire survey area, which results in large grid size and long computation time.

[0045] This application employs a local forward modeling strategy: for each undersampled emission source, a local grid is established and the electromagnetic field is solved only within a finite area around the undersampled emission source, rather than solving over the entire survey area grid. This local forward modeling strategy is based on the attenuation characteristics of the electromagnetic field in a conductive medium: during propagation, the energy of the electromagnetic wave is gradually absorbed by the medium, and the field value at a sufficiently far distance from the emission source has attenuated to a negligible level, making no substantial contribution to the response at the receiving point. A local grid is established for each undersampled emission source, and the three-dimensional electromagnetic response of that undersampled emission source is obtained through local solving, with independent calculations for each undersampled emission source. Since the local grid size is much smaller than the global grid, the computational cost of a single forward modeling operation is significantly reduced.

[0046] Centered on the current spatial location of the undersampled emission source, the local effective influence area is determined based on the skin depth. Skin depth is the distance an electromagnetic wave travels when its amplitude decays to approximately one-third of its original value while propagating in a conductive medium. Its value is related to the resistivity of the medium and the frequency of the electromagnetic wave: the higher the frequency, the smaller the skin depth, and the faster the electromagnetic field decays in space; the lower the resistivity, the smaller the skin depth, and the faster the decay. Therefore, for high-frequency components or low-resistivity media, a smaller local effective influence area can be selected; conversely, for low-frequency components or high-resistivity media, a larger area is required to meet accuracy requirements.

[0047] This embodiment employs a frequency-domain forward modeling method, setting multiple calculation frequencies covering the frequency band of a typical airborne electromagnetic system. The maximum skin depth corresponding to the lowest frequency is used as the primary basis for determining the local region's range, ensuring that the field value at the local region boundary has attenuated to a negligible level across all frequency components. Cutoff distances are determined along both the horizontal and vertical directions, forming a hexahedral local region centered on the transmitter location and bounded by the cutoff distance. The cutoff distance is taken as a preset multiple of the maximum skin depth corresponding to the transmitter location, ensuring that the electromagnetic field amplitude at the boundary is attenuated to below a preset threshold relative to the transmitter location.

[0048] After determining the spatial extent of the local effective influence region, a local mesh is established within this region. The cell size of the local mesh is consistent with the cell size of the central region of the global mesh to ensure the correct mapping of resistivity model parameters between different meshes. The local mesh covers the aforementioned hexahedral local region. The same cell size as the global mesh is used for the central region near the current emission source location to ensure computational accuracy. Near the boundary, the cell size can be appropriately increased to control the mesh size and reduce memory consumption and computational load.

[0049] In one example, the local grid is aligned with the central region of the global grid in both the horizontal and vertical directions, with a cell size of 60 meters, consistent with the basic cell of the global grid. The local region covers 240 meters horizontally and 300 meters vertically. Because the local grid's area is much smaller than the global grid's 6000m × 6000m × 2169m, its number of cells is much less than that of the global grid.

[0050] An extended element layer extends outward from the boundary of the local mesh, with the size of the extended element layer increasing by a preset multiple from the inside out. The purpose of the extended element layer is to absorb the electromagnetic wave energy propagating to the boundary without significantly increasing the computational load, attenuating the electromagnetic waves to a negligible level before reaching the boundary of the computational domain, thereby reducing the impact of boundary reflections on the accuracy of the internal solution. After setting the extended element layer, the electromagnetic field solution domain is extended outward from the physically relevant effective influence region by several element layers. When a truncation boundary condition is applied at the outer boundary, the influence of the boundary condition on the solution of the internal region is effectively suppressed because the electromagnetic waves have been sufficiently attenuated within the extended element layer.

[0051] Within the local grid, the total field is decomposed into a background field and a secondary field, which are calculated separately. The background field is the electromagnetic field generated by the current emission source under the background resistivity model (homogeneous half-space or layered medium, without three-dimensional anomalies), and can be efficiently calculated using analytical formulas or one-dimensional semi-analytical methods. The secondary field is the additional electromagnetic field generated by the difference between the actual resistivity model and the background resistivity model (i.e., the presence of three-dimensional anomalies).

[0052] The advantages of calculating the background field and the secondary field separately are as follows: the background field varies gently in space and can be calculated accurately and quickly using analytical methods; the secondary field is caused only by the anomalous body, varies gently over most of the local grid area, and only changes drastically near the anomalous body, resulting in a more stable numerical solution. Furthermore, the secondary field values ​​are very small at the local grid boundaries, making it easier to apply boundary conditions. In contrast, the total field exhibits singularities near the emission source, making it difficult to guarantee the accuracy of the numerical solution.

[0053] The background field and the solved secondary field are superimposed node by node to obtain the total electromagnetic response of the current undersampled transmitter under the current resistivity model, including the electric and magnetic field values ​​at each receiving point. The receiving point locations are distributed throughout the global measurement area, and the solution results of the local grid are mapped to each receiving point location through interpolation between grid nodes to obtain the response values ​​of each measurement time channel corresponding to the current transmitter.

[0054] In one example, the total electric field is separated from the background field and the secondary field using a method that distinguishes between the two fields. and total magnetic field Represented as: , , in, As the background electric field, As the background magnetic field, For secondary electric field, It is a secondary magnetic field; neglecting the displacement current, the diffusion equation of the secondary electric field is: , in, and The conductivity values ​​are those of the actual resistivity model and the background resistivity model, respectively. Permeability, For time.

[0055] For the above equation in a finite control volume Integrating and applying Gaussian divergence theorem, we obtain the governing equations for the finite volume: , in, For control body The boundary surface, Let be the unit outward normal vector on the boundary surface.

[0056] Each cell in the local mesh is treated as an independent control volume. A volume integral of the second-order electric field diffusion equation is performed on each control volume, and the curl term in the volume integral is transformed into a boundary surface integral using the Gaussian divergence theorem. After this transformation, the original partial differential equation becomes an algebraic conservation relation for the electromagnetic field quantities on each control volume. Repeating this process for all control volumes and assembling the equations yields a large-scale sparse linear system of equations, where the unknowns are the second-order electric field values ​​at each cell node. Solving this system using an iterative linear solver yields the second-order electric field distribution at each node of the local mesh, which is then combined with the background field to obtain the predicted response of the undersampled emission source.

[0057] Since subsequent sensitivity calculations require locating the parameters corresponding to the current local mesh in the global network model space, the index mapping relationship between local mesh elements and global mesh elements is recorded simultaneously with the local forward modeling. This mapping relationship is stored as an integer array, the length of which is the total number of elements in the current emitter's local mesh. The j-th element records the global element index in the global mesh corresponding to the j-th element of the local mesh. This index mapping is used in subsequent sensitivity calculations to quickly extract the local sensitivity row corresponding to the local mesh from the global network model space, avoiding out-of-bounds access and misaligned projection. The index mapping file for each undersampled emitter is named with its one-dimensional index in the complete emitter sequence for quick location and retrieval.

[0058] In step S103 of this application, the three-dimensional electromagnetic response is converted into an undersampled time-domain predicted response for each measurement time channel, and the complete time-domain predicted response corresponding to each emission source position is iteratively reconstructed in the sparse transform domain, including: The frequency domain response of each undersampled emission source is transformed into the frequency-time domain to obtain the undersampled time domain predicted response of each measurement time channel; For each measurement time channel, the undersampled time domain predicted response is arranged into a two-dimensional response matrix according to the spatial location of the emission source, the unsampled positions are set to zero, and a binary sampling position indication matrix of the same size is established. The curvelet transform is used as a sparse transform to iteratively reconstruct the two-dimensional response matrix for each measurement time channel. Each iteration includes: forward curvelet transform, thresholding the curvelet coefficients to zero according to the current threshold, and inverse curvelet transform. The threshold decreases with each iteration. The response value of the sampled position is forcibly replaced with the corresponding value in the original undersampled time domain predicted response according to the binary sampling position indication matrix. After iterative convergence, the reconstruction results of each measurement time channel are assembled in a predetermined order into the complete time-domain predicted response corresponding to each emission source position.

[0059] In some embodiments, the three-dimensional electromagnetic response of each undersampled emission source obtained through the local forward modeling calculation in step S102 is a frequency domain response, i.e., the electric and magnetic field values ​​obtained by solving Maxwell's equations at multiple preset calculation frequencies. However, the actual measurement of the time-domain airborne electromagnetic system is the secondary field attenuation response on each measurement time channel after the emission current is turned off, i.e., the time-domain signal. Therefore, it is necessary to convert the frequency domain response into the time-domain predicted response of each measurement time channel so that the time-domain predicted data and the observed data are directly comparable in the time domain.

[0060] The basic principle of frequency-time domain transformation is Fourier transform theory: the frequency domain response can be converted to the time domain through inverse Fourier transform. However, since the number of frequency points and time channels in actual calculations are finite discrete values, and the transmitted waveform of the airborne electromagnetic system is usually a step-off or pulse waveform, a transformation method suitable for step excitation is required. In this application, the embodiment uses cosine transform to achieve frequency-time domain transformation. For each undersampled transmitter, each receiver location, and each field component, its frequency domain response (including real and imaginary parts) at all calculated frequencies is subjected to cosine transform to obtain the undersampled time domain predicted response for each measurement time channel.

[0061] To suppress the Gibbs effect (oscillation artifact) caused by finite frequency bandwidth, window functions are applied to the low-frequency and high-frequency bands of the frequency domain response before transformation, so that the frequency domain response smoothly transitions to zero at the frequency band boundary. In this embodiment, the Hanning window function is used to weight the frequency domain response to reduce oscillations caused by the truncation effect.

[0062] After frequency-time domain transformation, the undersampled time-domain predicted responses of each undersampled transmitter are obtained for each measurement time channel. This means that only a portion of the transmitter locations (approximately 40%) have response values ​​across all measurement time channels; the response values ​​for the remaining transmitter locations have not yet been calculated. These undersampled time-domain predicted responses are then spatially reconstructed. Using a sparse transform domain iterative reconstruction method, the correlation of the electromagnetic response in spatial distribution is utilized to recover the complete time-domain predicted responses of all transmitter locations for each measurement time channel from the undersampled data.

[0063] For each measurement time channel, the undersampled time-domain predicted responses of each undersampled transmitter in that measurement time channel are arranged into a two-dimensional undersampled response matrix according to the spatial position of the transmitter in the two-dimensional plane of the survey area. The number of rows in this two-dimensional undersampled response matrix is ​​equal to the number of transmitter positions in one horizontal direction of the survey area, and the number of columns is equal to the number of transmitter positions in the other horizontal direction. The transmitter positions that have been calculated through forward modeling are filled with their true response values ​​in the two-dimensional undersampled response matrix, while transmitter positions that have not been calculated through forward modeling are filled with zero values. Through this rearrangement operation, the response of each measurement time channel is organized into a two-dimensional image, where the horizontal and vertical axes correspond to the two spatial dimensions of the survey area, and the pixel value corresponds to the response amplitude of the transmitter at that location in the current time channel.

[0064] Simultaneously, a binary sampling position indicator matrix with the exact same size as the two-dimensional response matrix is ​​established, where elements corresponding to sampled positions are set to 1, and elements corresponding to unsampled positions are set to 0. This binary sampling position indicator matrix remains unchanged throughout the reconstruction process and is used to mark which response values ​​are actually calculated and which need to be recovered through reconstruction.

[0065] Airborne electromagnetic responses exhibit significant spatial correlation: the changes in electromagnetic responses from adjacent emission source locations are continuous and gradual, without isolated abrupt changes. This spatial correlation makes the response compressible in the sparse transform domain, with its energy concentrated on a small number of large coefficients, while most coefficients are close to zero. This application employs curvelet transform as the sparse transform, which can more sparsely represent the electromagnetic response image with curved spatial characteristics, and is particularly suitable for the characteristic of curved singularities at the edges of anomalous bodies in airborne electromagnetic responses.

[0066] In the sparse transform domain, the energy of the complete response signal is concentrated in a small number of large coefficients, while the energy of artifacts caused by undersampling is dispersed among numerous small coefficients. By repeatedly performing an iterative process of forward transform, thresholding, inverse transform, and data backfilling, artifacts can be gradually eliminated while preserving the main energy of the signal, making the reconstruction result successively approximate the true complete time-domain predicted response.

[0067] At the start of each iteration, a curvelet forward transform is performed on the current two-dimensional response matrix to transform it to the curvelet domain, resulting in a curvelet coefficient matrix. Each coefficient in the curvelet coefficient matrix corresponds to a specific spatial scale and directional characteristic.

[0068] The current threshold is then determined based on the current iteration number. A larger threshold is initially set at the beginning of the iteration, retaining only the largest coefficient, corresponding to the main structure and spatial variation characteristics of the response. As the iteration progresses, the threshold is gradually decreased, restoring smaller coefficients to represent the local details and finer features of the response. Coefficients in the curvelet coefficient matrix whose absolute values ​​are below the current threshold are set to zero, while coefficients whose absolute values ​​are not lower than the current threshold remain unchanged, thus completing the thresholding process.

[0069] Then, the inverse curvelet transform is performed on the curvelet coefficient matrix after thresholding to transform the coefficients from the curvelet domain back to the spatial domain, resulting in the updated two-dimensional response matrix.

[0070] Finally, based on the binary sampling location indicator matrix, the elements in the updated two-dimensional response matrix corresponding to the sampled locations are forcibly replaced with the actual calculated values ​​at the corresponding locations in the original undersampled time-domain predicted response, while the reconstructed values ​​at unsampled locations remain unchanged. This forced backfilling operation ensures that the values ​​at the sampled locations are absolutely accurate throughout the entire iteration process, keeping the iteration always anchored to the true forward modeling results.

[0071] Repeat the above iterative process, decreasing the threshold each time, to quickly recover the main structure of the response in the early stages of the iteration and to finely recover the detailed features in the later stages. When the preset maximum number of iterations is reached or the relative change between two adjacent iterations is lower than the preset convergence threshold, the iteration is terminated, and the current two-dimensional response matrix is ​​the complete reconstruction result of the measurement time channel.

[0072] The spatial variation characteristics of the complete time-domain predicted responses of each measurement time channel differ: the responses of early time channels mainly reflect shallow information and exhibit gentle spatial variation; the responses of later time channels contain deep information and exhibit relatively complex spatial variation. Each measurement time channel independently performs the aforementioned iterative reconstruction, which can adapt to the spatial characteristics of different time channels. Initial thresholds and convergence criteria are set separately for each time channel to obtain optimal reconstruction accuracy.

[0073] After iterative convergence, a complete two-dimensional response matrix is ​​obtained for each measurement time channel, and all emission source locations have response values ​​in all time channels, thus realizing the recovery from undersampled data to the complete measurement area response.

[0074] In one embodiment, the three-dimensional electromagnetic response is converted into an undersampled time-domain predicted response for each measurement time channel. Let the first... The two-dimensional response matrix of each time channel is The set of undersampled locations is Ω, and the corresponding sampling operator is: Then, the data can be directly calculated to satisfy: , in, This is for predicting the response in the undersampled time domain.

[0075] Fill the results with zero values As the initial value, the h-th curve domain reconstruction iteration is expressed as: , in, and These are the forward curve transform and the inverse curve transform, respectively. For hard thresholding operators, the threshold is determined according to: , Decrease successively As the initial threshold, For the first The threshold of the next iteration Threshold decay factor, For the first The measurement time channel is at the first The two-dimensional response matrix obtained after iteration. For the first The measurement time channel is at the first The two-dimensional response matrix obtained after iteration. For the transpose of the sampling operator, For sampling operators, For unit operators, This is the curve coefficient matrix.

[0076] After steps S101 to S103, this application has been able to reconstruct the complete time-domain predicted response of all emission source locations from the forward modeling results of the undersampled emission sources. Based on this, the resistivity model needs to be continuously updated through inversion iterations to gradually approximate the observed data. The Gauss-Newton method is one of the most commonly used optimization methods for solving nonlinear inversion problems. Its core idea is to linearize the nonlinear forward modeling operator near the current resistivity model, transforming the nonlinear problem into a series of linear least-squares subproblems, and gradually approximating the true solution through iterative solutions.

[0077] This application first performs several iterations of full-order Gauss-Newton inversion. Full-order Gauss-Newton inversion refers to solving for model updates in the original high-dimensional model space (i.e., the parameter space composed of all grid cells), with each iteration constructing and solving the complete full-order Gauss-Newton equations. The purpose of full-order Gauss-Newton inversion is to ensure that the resistivity model converges sufficiently to the vicinity of the true solution, providing a high-quality source model for subsequent POD reduction inversion, i.e., the central model used to construct the POD snapshot. If the number of full-order iterations is insufficient, the model has not yet converged sufficiently, and the POD basis constructed based on this resistivity model cannot accurately reflect the structural characteristics of the true model, and the reduction inversion may deviate from the correct solution; if the number of full-order iterations is too large, the computational advantages of full-order Gauss-Newton inversion are not fully utilized, reducing overall computational efficiency.

[0078] In the Kth full-order iteration, the current resistivity model is first read (in the first iteration, the current resistivity model is the prior resistivity model). The complete time-domain predicted response under the current resistivity model is then invoked. Specifically, based on the determined set of undersampled emission sources, local three-dimensional forward modeling is performed on each undersampled emission source to obtain the frequency-domain electromagnetic response of each undersampled emission source; then, through frequency-time domain transformation and sparse domain iterative reconstruction, the complete time-domain predicted response vectors of all emission source positions across all measurement time channels are recovered.

[0079] The full-order weighted residuals are calculated based on the difference between the observed data and the predicted response, combined with the data covariance matrix. The full-order residual vector is defined as the difference between the observed data and the current resistivity model's predicted response, with each element corresponding to the difference between the observed and predicted data. The larger the residual, the further the predicted response of the current resistivity model deviates from the measured data, and the greater the correction amount required for the resistivity model.

[0080] The diagonal elements of the data covariance matrix represent the estimated error variance of each observation. The residuals are weighted using the square root inverse of the data covariance matrix to obtain the full-order weighted residuals. In the full-order weighted residuals, data with high quality and low noise receive greater weights and play a more important role in the inversion; data with low quality and high noise receive less weights, thus reducing their impact on the inversion results. This weighting method makes the inversion focus more on reliable data, improving the robustness of the inversion.

[0081] The full-order weighted sensitivity matrix (also known as the Jacobian matrix) is one of the most crucial mathematical objects in the inversion process. Its elements are defined as the partial derivatives of each predicted response with respect to each resistivity model parameter, representing the change in the i-th predicted data when the j-th resistivity model parameter undergoes a small change. The full-order weighted sensitivity matrix establishes a linearized mapping relationship between the resistivity model space and the data space, determining the computational accuracy of the resistivity model update direction.

[0082] In this application, the full-order weighted sensitivity matrix is ​​calculated using the adjoint forward modeling method. This method can obtain the sensitivity of an observation data point to all resistivity model parameters through a single adjoint forward modeling, and its computational efficiency is much higher than that of the element-by-element perturbation method. Specifically, for each emission source, the adjoint field equation is solved once, and the time-domain convolution result of the adjoint field and the forward modeling field is mapped to each grid cell to obtain the sensitivity values ​​of all measured time-channel responses of the emission source to all resistivity model parameters. It should be noted that the full-order weighted sensitivity matrix is ​​composed of multiple full-order weighted sensitivities arranged together.

[0083] Define full-order weighted residuals and full-order weighted sensitivity matrix : , , in, This is the unweighted sensitivity matrix.

[0084] Because this application employs a local forward modeling strategy, the sensitivity of each emission source is also stored in a local window format according to the emission source location. The number of rows in the local weighted sensitivity of each emission source is equal to the number of time channels of that emission source, and the number of columns is equal to the number of cells in the local grid of that emission source. When assembling the full-order weighted sensitivity matrix, the local-to-global index mapping relationship needs to be used to extract the corresponding rows from each local weighted sensitivity matrix and place them in the correct positions of the full-order weighted sensitivity, ensuring that the row order of the matrix composed of full-order weighted sensitivity is strictly consistent with the arrangement order of the predicted response vector.

[0085] After obtaining the full-order weighted residuals and the full-order weighted sensitivity matrix, the full-order Gauss-Newton equation can be constructed, and the direction of resistivity model update can be obtained by solving the full-order Gauss-Newton equation.

[0086] The gradient of the objective function is determined by the full-order weighted residual and the full-order weighted sensitivity matrix. Specifically, it is the transpose of the full-order weighted sensitivity matrix multiplied by the full-order weighted residual, which reflects the direction of the fastest rise of the objective function in the resistivity model space. Its negative direction is the direction of the fastest fall.

[0087] The Gauss-Newton approximation Hessian matrix consists of two parts: the transpose of the full-order weighted sensitivity matrix multiplied by itself, plus the product of the regularization parameter and the inverse of the model covariance matrix. The transpose of the full-order weighted sensitivity matrix multiplied by itself describes the local curvature information of the objective function, ensuring that the update direction includes not only gradient information but also curvature correction, resulting in a faster convergence speed than the steepest descent method alone. The regularization term improves the ill-conditioned conditions of the Gauss-Newton approximation Hessian matrix, ensuring the numerical stability of the linear equation system.

[0088] gradient of the objective function and the Gauss-Newton approximation Hessian matrix They are represented as follows: , , After obtaining the full-order weighted residuals and the full-order weighted sensitivity matrix, solving the full-order Gauss-Newton model yields the resistivity model update direction. The full-order Gauss-Newton equation is expressed as: , in, This provides a direction for updating the resistivity model.

[0089] The full-order Gauss-Newton equations are a system of linear equations with the Gauss-Newton approximation Hessian matrix as the coefficient matrix and the gradient of the objective function as the right-hand side. Solving the linear equations yields the direction of resistivity model updates, which indicates how the resistivity values ​​of each grid cell should be modified based on the current resistivity model to reduce the objective function.

[0090] After obtaining the direction of resistivity model updates, the current resistivity model cannot be directly modified by adding the update direction to create the new resistivity model, as the update step size might be too large, causing the objective function to increase. Therefore, a line search is needed to determine the optimal step size to ensure a sufficient decrease in the objective function.

[0091] The specific process of the full-order line search is as follows: Using a step size of 1 as the initial step size, a full-order trial resistivity model is generated along the resistivity model update direction. The complete time-domain predicted response of this full-order trial resistivity model is calculated, and its objective function value is evaluated. It is then determined whether the objective function value of this full-order trial resistivity model satisfies the Armijo sufficient descent condition: the new objective function value is not greater than the sum of the current objective function value and the preset proportion of the step size multiplied by the gradient direction projection. Armijo refers to the fact that the objective quasi-descent of the new trial resistivity model must not be less than the theoretical maximum descent of the preset proportion.

[0092] If the full-order trial resistivity model meets the descent condition, it is accepted as the updated resistivity model for this iteration; if it does not meet the condition, the full-order trial resistivity model is regenerated after reducing the step size by a preset ratio, the predicted response and objective function value are recalculated and re-evaluated, and the process is repeated until the descent condition is met or the maximum number of attempts is reached.

[0093] The first full-order trial resistivity model that satisfies the descent condition is marked as the accepted full-order trial resistivity model, and the current full-order trial resistivity model, predicted response, and objective function value are updated. This acceptance marking is crucial; only resistivity models that satisfy the descent condition have physical rationality and numerical stability. The source resistivity model selected in subsequent POD descent stages must be an accepted full-order trial resistivity model that has been verified through line search.

[0094] After completing one full-order iteration, it is determined whether the current inversion state meets the preset order reduction switching conditions. The switching conditions include, but are not limited to: reaching the preset number of full-order iterations (2 in this embodiment), the relative change of the objective function being lower than a preset threshold, and the normalized root mean square error being lower than a preset threshold. When any condition is met, the full-order inversion is terminated, and the POD order reduction stage begins.

[0095] At this point, the last full-order trial resistivity model accepted by the line search in the full-order stage is selected as the POD source resistivity model. This POD source resistivity model is located in the original high-dimensional resistivity model space and represents the current best resistivity distribution estimate after the full-order inversion convergence. The reason for selecting the last accepted updated resistivity model rather than the latest calculated but unaccepted trial resistivity model is that only the accepted full-order trial resistivity model satisfies the objective function descent condition.

[0096] After the full-order Gaussian-Newton inversion in step S104, the resistivity model has initially converged to near the true solution. If inversion iterations continue in the full-order space, each iteration still requires constructing and solving a system of tens of thousands of linear equations, resulting in a considerable computational burden. However, in the later stages of the inversion iteration, the update of the resistivity model mainly involves fine-tuning of the local structure and refining of details, while the overall structure of the resistivity model tends to stabilize. At this point, the update of the resistivity model is actually constrained to a low-dimensional subspace much smaller than the dimension of the full-order space.

[0097] This application employs an intrinsic orthogonal decomposition method to reduce the order of the model space. The core idea is to extract an orthogonal basis that can characterize the main changes in the model with the fewest possible modes by statistically analyzing several possible model states (i.e., snapshots). This orthogonal basis is then used as the basis of a low-dimensional subspace, projecting the resistivity model update problem from the original high-dimensional space onto this low-dimensional space for solution, thereby significantly reducing the solution scale of the linear equation system.

[0098] When the order reduction switching condition in step S104 is met, the full-order inversion terminates. At this time, the last full-order trial resistivity model accepted by the line search in the full-order stage is read from the storage medium and used as the POD source resistivity model.

[0099] The POD source resistivity model is the central point for subsequent snapshot construction, and its quality directly determines the ability of the reduced-order basis to characterize the features of the real resistivity model. If the POD source resistivity model already reflects the real subsurface structure well, then the snapshot constructed based on it can effectively cover the possible range of variations in the real resistivity model.

[0100] Identify anomalous regions in the POD source resistivity model that deviate from the background resistivity model. Specifically, calculate the logarithmic deviation of the resistivity of each grid cell relative to the background resistivity (i.e., the logarithmic value of the resistivity of each grid cell minus the logarithmic value of the background resistivity). Cells with an absolute value of this logarithmic deviation greater than a preset threshold are marked as candidate cells. Then, perform spatial connectivity analysis on the candidate cells, and mark the set of spatially connected candidate cells as an anomalous region.

[0101] When no grid cell meets the preset threshold, it is determined that no obvious abnormal region has been detected. In this case, the entire grid model region is used as the candidate region to ensure that the subsequent snapshot construction will not fail due to the absence of abnormal regions.

[0102] Multiple snapshots are constructed centered on the POD source resistivity model. Each snapshot is a resistivity model vector obtained by applying a specific type of perturbation to the POD source resistivity model. The purpose of snapshot construction is to generate a set of samples that can cover various resistivity variation patterns that may appear in the inversion model, so that the subsequently extracted POD basis can effectively characterize the low-dimensional structure of the model space.

[0103] This application embodiment designs three perturbation modes: (1) Background amplitude perturbation: A random perturbation following a normal distribution is applied to all grid cells (i.e., background regions) located outside the anomalous region in the POD source resistivity model, causing their resistivity values ​​to change randomly within a certain range. The standard deviation of the background amplitude perturbation is taken as the first preset value, so that the background region produces small fluctuations, which is used to enhance the ability of the POD basis to express changes in the background field.

[0104] (2) Abnormal amplitude perturbation: A random perturbation following a normal distribution is applied to all grid cells located within the abnormal region in the POD source resistivity model, causing the resistivity value of the abnormal region to change randomly within a certain range. The standard deviation of the abnormal amplitude perturbation is taken as a second preset value, which is greater than the standard deviation of the background amplitude perturbation, so as to reflect the physical characteristic of the large range of resistivity changes in the abnormal body.

[0105] (3) Anomaly position perturbation: The anomaly region is translated by one or more mesh elements along the three coordinate directions of the 3D mesh to generate a snapshot of the anomaly's position shift. For the portion that extends beyond the mesh boundary after translation, it is truncated to the mesh boundary. This perturbation method is used to simulate the uncertainty of the spatial position of the anomaly and enhance the ability of the POD basis to characterize the changes in the anomaly's position.

[0106] Each model snapshot can contain only one of the three perturbation modes mentioned above, or it can contain a combination of multiple modes. By combining different perturbation types and perturbation parameters, a diverse set of snapshots can be generated.

[0107] After constructing all snapshots, each snapshot is arranged into a one-dimensional column vector, the length of which is equal to the total number of parameters in the regular mesh model. All snapshot vectors are then arranged column-wise to form a snapshot matrix. Each column of the snapshot matrix corresponds to a snapshot, and each row corresponds to the resistivity value of a mesh cell in different snapshots.

[0108] Calculate the mean vector for all snapshots, i.e., take the average of each row of the snapshot matrix to obtain a mean vector. Subtract the mean vector from each snapshot vector to obtain the mean-reduced snapshot matrix. The mean-reduction process makes the mean of each row of the snapshot matrix zero, eliminating the common background part in the snapshots, and allowing subsequent singular value decomposition to focus more on the differences and variation characteristics of the resistivity model.

[0109] The mean snapshot and the mean-reduced snapshot matrices are as follows: , , in, This is a snapshot of the mean. This is a snapshot matrix with the mean removed. For the total number of snapshots, This is the s-th snapshot vector. For the first A snapshot vector.

[0110] The mean-free snapshot matrix is ​​decomposed into singular value decomposition (SVD), resulting in a product of three matrices: a left singular vector matrix, a singular value diagonal matrix, and a right singular vector matrix. Each column of the left singular vector matrix is ​​a left singular vector, the length of which is equal to the total number of parameters in the regular grid model. The diagonal elements of the singular value diagonal matrix are singular values, arranged in descending order. The square of each singular value represents the energy share of the corresponding mode in the snapshot matrix.

[0111] Snapshot matrix with mean removed Singular value decomposition is represented as: , in, It is a left singular vector matrix. For singular value diagonal matrices and It is a right singular vector matrix.

[0112] The magnitude of singular values ​​directly reflects the importance of the corresponding mode: the larger the singular value, the more dominant the resistivity mode is in the snapshot set and the greater its contribution to characterizing the changes in the resistivity model; the smaller the singular value, the smaller the contribution of the resistivity mode to the changes in the resistivity model, and it can be ignored without significantly affecting the reconstruction accuracy.

[0113] The proportion of the sum of squares of the first r singular values ​​to the sum of squares of all singular values ​​is calculated as the cumulative energy percentage of the first r modes. This proportion reflects the percentage of the total energy of the snapshot matrix that can be preserved when using the first r left singular vectors as the basis. The smallest integer r that satisfies the requirement that the cumulative energy percentage is not less than a preset energy preservation threshold is selected as the reduced dimension.

[0114] The first r columns of the left singular vector matrix are taken as the POD basis in the internal cell space of the regular grid. The column vectors of this POD basis span an r-dimensional subspace in the internal cell space of the regular grid. The internal cell space of the regular grid refers to the high-dimensional vector space formed by defining the resistivity model parameters on a regular three-dimensional grid covering the entire survey area. Subsequent resistivity model updates will be performed in this r-dimensional subspace.

[0115] Choose the minimum reduced dimension r based on the cumulative energy of the squared singular values: , in, To maintain the energy threshold, the first r columns of left singular vectors are taken to form the POD basis. It is the kth singular value.

[0116] The obtained POD basis is constructed on the internal cell space of a regular mesh, with the number of rows equal to the total number of parameters in the regular mesh model and the number of columns equal to the reduced dimension. However, the actual solution for the full-order Gaussian-Newton inversion is performed in the internal inversion parameter space, which has a model covariance transformation relationship with the internal cell space. Simultaneously, the calculated weighted sensitivity matrix is ​​also defined in the internal inversion parameter space. If the POD basis in the internal cell space of the regular mesh is directly used for the reduced-order projection, spatial inconsistency problems will arise: mesh dimension mismatch, inconsistent parameter order, and inconsistent physical meanings. This leads to the reduced objective function being inequivalent to the original full-order objective function, and the inversion iteration may diverge or converge to an incorrect solution.

[0117] Therefore, a two-level mapping is required to transform the POD basis from the internal cell space of the regular grid to the internal inversion parameter space. The first-level mapping involves filtering and rearranging the regular grid cells into the internal inversion cells, and the second-level mapping involves transforming the internal cell space into the internal inversion parameter space.

[0118] The regular grid is a complete three-dimensional grid covering the entire survey area, including all cells in the air layer, background region, and subsurface anomaly region. The regular grid contains regular grid cells. However, during the inversion process, typically only cells within a certain depth range that may contain anomalies are updated; cells in the air layer and at the bottom of the deep half-space are not involved in the inversion. These cells involved in the inversion are called internal inversion cells, and their number is less than or equal to the total number of regular grid cells.

[0119] The first-level mapping is achieved through a constraint mapping matrix. Based on the coordinate correspondence between regular grid cells and internal inversion cells, a constraint mapping is established from the regular grid index to the internal inversion cell index, implemented in one example using a correspondence table. The correspondence table records the original position of each internal inversion cell in the regular grid and its sequence number in the internal inversion cell array. Using this correspondence table, the POD base of the regular grid is row-filtered, retaining rows corresponding to internal inversion cells and discarding rows corresponding to non-inversion cells (such as air layers), and the rows are rearranged according to the unified storage order of the internal inversion cells.

[0120] After the first-level mapping, the POD basis in the internal cell space is obtained. Consistency checks are required before and after the mapping: verifying whether the number of rows after mapping equals the total number of internal inversion cells, whether each element is a finite value, and whether the index range is within the legal boundaries. If any condition is not met, stop entering the reduced-order iteration and output diagnostic information or revert to the full-order iteration.

[0121] After the first-level mapping, the POD basis still resides in the inner cell space. The actual update of the full-order Gaussian-Newton inversion takes place in the inner inversion parameter space, and the two are linked by the square root of the model covariance operator. Therefore, the mapping equation between the POD basis in the inner cell space and the reduced-order basis in the inner inversion parameter space is: , For the POD basis in the internal unit space, This serves as a reduced-order basis for the internal inversion parameter space. It should be noted that the internal cell space refers to the physical space spanned by defining resistivity values ​​only on the grid cells involved in the inversion, after removing non-inversion cells (such as air layers and the bottom of deep half-spaces) from the regular grid covering the entire survey area. The internal inversion parameter space refers to the unconstrained mathematical parameter space obtained by applying a model covariance transformation to the resistivity physical space. The resistivity physical space is a generalized space spanned by the resistivity of all grid cells (including non-inversion cells such as air layers).

[0122] After the second-level mapping, the reduced basis has been completely transformed from the internal cell space to the internal inversion parameter space, resulting in the reduced basis of the internal inversion parameter space.

[0123] The consistency of the two-space mapping is verified by the following formula: , in, For consistency error, when If the threshold is exceeded, or if a non-numerical value or mapping dimension is inconsistent, the process will not proceed to the reduced-order iteration, and diagnostic information will be output or the process will revert to the full-order iteration.

[0124] After completing the two-level mapping, a reduced-order basis is obtained in the internal inversion parameter space. This reduced-order basis spans an r-dimensional low-dimensional subspace, and subsequent model updates will be performed in this low-dimensional subspace.

[0125] The order reduction inversion is: The model update increment of dimension is constrained to be a linear combination of the columns of the reduced-order basis, i.e. ,in Update the r-dimensional reduction coefficient vector. For reducing order basis, This represents the update amount for the reduced-order model. Thus, what was originally required to solve... The problem of solving a system of r-dimensional linear equations is transformed into solving only a system of r-dimensional linear equations. Because much smaller (In this embodiment) Typically, it is tens of dimensions, while (For systems of linear equations with tens of thousands of dimensions), the solution size is significantly reduced, and computational efficiency is significantly improved.

[0126] Using the reduced-order basis, the full-order weighted residuals and full-order weighted sensitivity matrices are projected onto the reduced-order space, and the reduced-order Gauss-Newton equations are solved to obtain the update direction of the reduced-order model.

[0127] The process involves projecting the full-order weighted sensitivity matrix onto the reduced-order space. The number of rows in the full-order weighted sensitivity matrix equals the total number of observed data points, and the number of columns equals the total number of internal inversion units. Multiplying the full-order weighted sensitivity matrix by the reduced-order basis yields the reduced-order weighted sensitivity matrix, whose number of rows still equals the total number of observed data points, and whose number of columns equals the reduced-order dimension. Each column of the reduced-order weighted sensitivity matrix corresponds to the degree of influence of a POD mode (multiple POD modes arranged in columns form a POD basis) on the predicted data. The reduced-order weighted sensitivity matrix describes the impact of changes in the low-dimensional reduced-order coefficients on the predicted response and is the core object for constructing the reduced-order Gauss-Newton equation.

[0128] When the full-order weighted sensitivity is stored in a local window format according to the emission source location, it needs to be processed according to the index mapping from the local grid cells to the global grid cells corresponding to each emission source. Specifically, for each observation data point, the corresponding local row is extracted from the full-order weighted sensitivity matrix using the index mapping from the local grid cells to the global grid cells corresponding to the emission source. This local row is then multiplied by the corresponding row in the reduced-order basis, constructing the reduced-order weighted sensitivity matrix row by row. This processing method avoids the out-of-bounds and misalignment problems caused by directly using the local sensitivity array as a global vector for calculation.

[0129] Multiplying the reduced-order weighted sensitivity matrix by the reduced-order weighted residual yields the reduced-order negative gradient vector. The left-hand coefficient matrix of the reduced-order Gauss-Newton equation is formed by multiplying the transpose of the reduced-order weighted sensitivity matrix with the product of the regularization parameter and the transpose of the reduced-order basis with the product of the reduced-order basis. The reduced-order Gauss-Newton equation is constructed by using the reduced-order negative gradient vector as the right-hand term.

[0130] Reduced right-hand term And the reduced-order Gaussian-Newton approximation Hessian matrix They are respectively: , , in, For the right-hand side of the reduced order term, For the reduced-order Gaussian-Newton approximation Hessian matrix, For the transpose of a reduced basis, It is the full-order negative gradient vector. This is a reduced-order weighted sensitivity matrix.

[0131] In the reduced-order space, the coefficient matrix of the reduced-order Gauss-Newton equation consists of two parts: The first part is the product of the transpose of the reduced-order weighted sensitivity matrix and the reduced-order weighted sensitivity matrix. This part describes the local curvature information of the objective function in the reduced-order space, whose dimension is the product of the reduced-order dimension and the reduced-order dimension.

[0132] The second part is the regularization term, which is the regularization parameter multiplied by the product of the transpose of the reduced-order basis and the reduced-order basis. Since the reduced-order basis does not satisfy the orthogonality condition in the internal inversion parameter space, the regularization term must be explicitly calculated as the transpose of the reduced-order basis multiplied by itself and then by the regularization parameter, rather than simplifying to the regularization parameter multiplied by the identity matrix. This regularization term maintains the formal consistency of the model constraint terms in the original objective function, ensuring that the reduced-order inversion maintains the same mathematical form as the original full-order objective function.

[0133] Adding the two parts together forms the coefficient matrix of the reduced-order Gauss-Newton equation. The right-hand side of the reduced-order Gauss-Newton equation is the reduced-order right-hand side. This equation is a system of linear equations with dimension equal to the product of the reduced-order dimension and the reduced-order dimension.

[0134] The reduced-order Gauss-Newton equation is expressed as: , in, This is the vector of reduced-order coefficients. To reduce the update volume, For weighted residuals, The coefficient matrix, This is a reduced-order weighted sensitivity matrix. This is the transpose of the reduced-order weighted sensitivity matrix. This is the current internal inversion parameter vector. This represents the right-hand side of the reduced-order linear system. The reduced-order linear system is solved using the conjugate gradient method with a maximum of 2000 iterations, ensuring that the squared relative residual is no greater than 10⁻⁻⁶. 7 As a stopping condition.

[0135] The reduced-order Gauss-Newton equation has a much smaller dimension than the full-order Gauss-Newton equation, significantly reducing the solution complexity. Solving this equation using the conjugate gradient method or other symmetric positive definite linear solvers yields the update direction of the reduced-order coefficients. The dimension of the update direction is equal to the reduced-order dimension, and its physical meaning is: based on the current reduced-order coefficients, how should the coefficients of each POD basis be adjusted to reduce the objective function?

[0136] Before solving, a consistency check is required: verify that the number of rows in the reduced-order weighted sensitivity matrix equals the total number of observations, and that the number of columns equals the reduced-order dimension; verify that the dimension of the reduced-order coefficient matrix is ​​the reduced-order dimension multiplied by the reduced-order dimension. If any condition is not met, stop the reduced-order iteration and revert to the full-order iteration.

[0137] After the solution is obtained, the subsequent reduction line search will generate trial reduction coefficients along the reduction coefficient update direction, and then convert them into a reduced trial resistivity model through the reduction basis, calculate the predicted response and evaluate the objective function value.

[0138] In step S108, the order-reduction model update direction is used to generate a reduced-order trial resistivity model through order-reduction line search, and the complete time-domain predicted response of the reduced-order trial resistivity model is calculated; only reduced-order trial resistivity models that satisfy the objective function value descent condition are accepted as the current inversion model, including: With a step size of 1 as the initial value, a trial reduction coefficient is generated along the update direction of the reduction model. Based on the linear combination of the current internal inversion model, the reduction basis, and the trial reduction coefficient, a trial internal inversion model is generated. Based on the mapping relationship between the internal inversion parameters and the resistivity model parameters, the trial internal inversion model is converted into a reduced-order trial resistivity model. Calculate the complete time-domain predicted response of the reduced-order trial resistivity model; The objective function value is calculated based on the difference between the observed data and the predicted response in the complete time domain, and it is determined whether the objective function value satisfies the descent condition relative to the objective function value of the current inversion model. If the descent condition is met, the trial resistivity model is accepted as the current inversion model, and the corresponding full time domain prediction response, objective function value, and normalized root mean square error are updated. If the descent condition is not met, the trial order reduction coefficient, trial internal inversion model, and reduced-order trial resistivity model are regenerated after reducing the step size by a preset ratio. The full time domain prediction response is recalculated and re-evaluated until the descent condition is met or the maximum number of attempts is reached. If no reduced-order resistivity model satisfies the descent condition after reaching the maximum number of attempts, the current inversion model remains unchanged, and the current line search is terminated.

[0139] The basic idea of ​​the order reduction line search is to generate a series of candidate models along the update direction of the order reduction model, starting from the initial step size, calculate the complete predicted response and objective function value of each candidate model, until a candidate model that meets the descent condition is found, and accept it as the update result of this iteration.

[0140] The essential difference between reduced-order line search and full-order line search lies in the fact that full-order line search directly adjusts the resistivity values ​​of each grid cell in the high-dimensional model space, while reduced-order line search adjusts the coefficients of each POD mode in the low-dimensional reduced-order space, and then maps them to the high-dimensional model space through the reduced-order basis to obtain the candidate resistivity model. The former modifies tens of thousands of independent parameters, while the latter only adjusts dozens of reduced-order coefficients.

[0141] The order reduction line search starts with a step size of 1 and generates trial order reduction coefficients along the obtained order reduction model update direction. Specifically, the current order reduction coefficient is added to the product of the step size and the order reduction model update direction to obtain the trial order reduction coefficients.

[0142] In each descent line search, the step size decreases sequentially according to a preset sequence. In this embodiment, the step size sequence is 1, 0.5, 0.25, and 0.125, meaning that the step size for each attempt is half of the previous one. Corresponding trial descent coefficients are generated using each step size value, corresponding candidate models are generated, and the objective function value is calculated until a step size that satisfies the descent condition is found or the maximum number of attempts is reached.

[0143] For each set of trial reduction coefficients, they need to be converted into a resistivity model that can be used in three-dimensional forward modeling. The conversion process consists of two steps: The first step is to multiply the trial reduction coefficients by the reduction basis to obtain the update amount of the trial internal inversion parameters. This update amount is then added to the current internal inversion parameters to obtain the trial internal inversion parameters.

[0144] The second step involves transforming the experimental internal inversion parameters into a reduced-order experimental resistivity model based on the model transformation relationship between the established resistivity model parameters and the internal inversion parameters. Specifically, this involves adding the product of the prior resistivity model and the square root of the model covariance operator with the experimental internal inversion parameters. After this transformation, a resistivity value distribution covering the entire three-dimensional grid of the survey area is obtained, which can then be used for forward modeling calculations.

[0145] After generating the reduced-order experimental resistivity model, its physical rationality needs to be checked. This mainly includes: whether the resistivity values ​​of each grid cell are all positive; whether the resistivity values ​​are within a preset reasonable range (e.g., between 0.1 Ω·m and 10000 Ω·m); and whether there are any non-numerical values. If any of the above checks fail, the reduced-order experimental resistivity model is rejected directly, and the next step is attempted.

[0146] For each reduced-order trial resistivity model, its complete time-domain predicted response is calculated. This calculation process is exactly the same as the full-order inversion stage: based on the determined set of undersampled emitters and their index mapping, a local three-dimensional forward model is performed on each undersampled emitter to obtain the frequency domain response of each undersampled emitter; then, through frequency-time domain transformation and curvelet domain iterative reconstruction, the complete time-domain predicted response of all emitter locations on all time channels is recovered.

[0147] The dimensions and order of the complete time-domain predicted response are completely consistent with the observed data, and can be directly used for evaluating the objective function value. The reason for re-executing the complete forward modeling and response reconstruction process for each order reduction line search is to ensure the accuracy of the objective function evaluation—avoiding the use of linearized approximate predicted values ​​to replace the true forward modeling results, and preventing the acceptance of unreliable models due to approximation errors.

[0148] After obtaining the complete time-domain predicted response of the reduced-order test resistivity model, the objective function value of the model is calculated. The objective function consists of two parts: a data-weighted fitting term and a model constraint term. The data-weighted fitting term is the sum of squares of the differences between the observed data and the predicted response, weighted by the data covariance. The model constraint term is the product of the square of the L2 norm of the test internal inversion parameter vector and the regularization parameter.

[0149] Determine whether the reduced-order trial resistivity model satisfies the sufficient descent condition. This condition means that the objective function value of the reduced-order trial resistivity model decreases significantly relative to the objective function value of the current resistivity model, and the actual decrease is not less than a preset proportion of the decrease predicted by the linear approximation.

[0150] If the sufficient descent condition is met, the trial resistivity model is accepted, and the search for the descent line is terminated. If not, the step size is reduced by a preset ratio, the trial descent coefficients, the trial internal inversion model, and the descent trial resistivity model are regenerated, the predicted response and objective function values ​​are recalculated, and the descent condition is re-evaluated. This process is repeated until a step size that satisfies the descent condition is found or the maximum number of attempts is reached.

[0151] When a reduced-order trial resistivity model satisfies the sufficient descent condition, it is accepted as the current inversion model, and the following state update is performed: Update the order reduction coefficients before order reduction to the order reduction coefficients corresponding to the experimental resistivity model.

[0152] Update the current resistivity model to this reduced-order trial resistivity model.

[0153] Update the current predicted response to the full time-domain predicted response of the experimental resistivity model.

[0154] Update the current objective function value to the objective function value of this reduced-order trial resistivity model.

[0155] Update the normalized root mean square error (RMSE). RMSE is a commonly used metric for measuring the goodness of fit of data. It is calculated by summing the squares of the differences between the observed data and the predicted response, dividing by the standard deviation of each data point, then dividing by the total number of valid data points and taking the square root. A smaller RMSE indicates that the predicted response of the current model is closer to the observed data. When the RMSE decreases below the preset target value, it can be used as one of the conditions for terminating the inversion.

[0156] If no reduced-order trial resistivity model that meets the descent condition is found within the maximum number of attempts, the model will not be updated in this iteration. The current model will remain unchanged and the line search will be terminated. Then, it will be determined whether the termination condition is met. If not, the POD basis will be reconstructed or the iteration will be rolled back to the full-order iteration.

[0157] In one example, in the k-th reduction iteration, a reduction line search is performed based on the update direction of the reduction model. Indicates the first The trial step size for the second-order reduction line search is used to construct the trial reduction coefficient and the reduced-order trial resistivity model according to the following formula. When the... When the trial satisfies the descent condition of the objective function, the corresponding trial result is used as the model for the (k+1)th iteration.

[0158] , , in, For the first In the second reduction iteration, the first The second attempt at order reduction line search is to test the order reduction coefficient. For the first The current reduction coefficient in the next reduction iteration For the first The direction of updating the order reduction model in the next order reduction iteration. For the first In the second reduction iteration, the first A reduced-order tentative resistivity model attempted by a secondary line search. This is the currently accepted reduced-order trial resistivity model.

[0159] For each reduced-order trial resistivity model, steps S102 to S103 are called back to calculate and reconstruct the complete time-domain predicted response at the same undersampled emission source location, and the original full-order objective function is re-evaluated; it is accepted only if the descent condition is met, otherwise the step size s is reduced.

[0160] After accepting the time-domain trial resistivity model, the predicted response, objective function, and normalized root mean square error are updated. Normalized root mean square error Defined as: , in, For the first The standard deviation of each observation data point The total number of valid data. For the first The measured values ​​of each observation data point For the first The predicted response value corresponds to each observation data point. If the termination condition is not met, return to step S106; if the termination condition is met, output the inversion model. The termination condition includes at least one of the following: the normalized root mean square error reaches the target value, the relative change of the objective function is lower than a threshold, the reduced negative gradient norm is lower than a threshold, the update amount is lower than a threshold, and the maximum number of iterations is reached; if the POD reconstruction error, basis mapping error, or reduced linear solution residual does not meet the threshold, perform a full-order backtracking or reconstruct the POD basis.

[0161] This application uses a three-dimensional time-domain airborne electromagnetic model with a low-resistivity cube set in a uniform half-space as an example to verify the POD order reduction inversion method. Figure 2 This is a schematic diagram illustrating the co-mapping between the POD basis of the regular mesh and the internal inversion parameter space. First, the POD basis is constructed in the resistivity space of the regular mesh based on a model snapshot. Figure 2 (a) represents the POD basis of the regular mesh; subsequently, based on the spatial correspondence between the regular mesh elements and the internal inversion elements, the POD basis of the regular mesh is subjected to effective element selection and parameter order rearrangement, so that its row count is converted from the number of regular mesh elements to the number of internal inversion elements, see... Figure 2 In (b), not all regular mesh elements participate in the inversion; only valid elements are extracted and rearranged according to the order of their internal inversion parameters. See [reference needed]. Figure 2 (c) represents the POD basis in the rearranged internal cell space. Based on this, and combining the parameter transformation corresponding to the square root of the model covariance, a reduced-order basis is formed that is consistent with the internal inversion parameter dimension, cell order, and parameterization method, as shown in [reference needed]. Figure 2 In the diagram, (d) represents the reduced-order basis of the internal inversion parameter space. Finally, the full-order weighted sensitivity matrix is ​​projected onto the reduced-order basis to obtain a reduced-order weighted sensitivity matrix with the number of columns equal to the reduced-order dimension r, as shown in the diagram. Figure 2 In this context, (e) represents the full-order weighted sensitivity matrix. By using a reduced-order basis, the full-order weighted sensitivity matrix is ​​transformed into... Figure 2 In the equation (f), the reduced-order weighted sensitivity matrix is ​​represented, and the updated direction of the reduced-order model is solved in the r-dimensional reduced-order space. Figure 2 The shaded cells in the diagram represent the corresponding mesh components that are retained and participate in model updates during the mapping process. This mapping avoids the inconsistency between the POD basis of the regular mesh and the internal inversion parameter space in terms of model dimensions and cell arrangement order.

[0162] Figure 3 The actual resistivity model used in this embodiment is set to a background resistivity of 100 Ω·m and a low-resistivity anomaly resistivity of 1 Ω·m; the model has a side length of approximately 240 m, is located in the center of the survey area, and has a top burial depth of approximately 90 m. First, the time-domain airborne electromagnetic response is calculated based on the actual resistivity model, and a Gauss-Newton inversion is performed using a uniform half-space model as the initial inversion model. The initial RMS is 2.96. After completing the second full-order iteration, a snapshot is constructed using the accepted full-order trial resistivity model as the POD source resistivity model. The snapshot matrix is ​​then subjected to mean removal and singular value decomposition, and the POD basis is selected according to the energy preservation criterion, thus entering the order reduction inversion stage. After subsequent order reduction iterations, the RMS (normalized root mean square error) is reduced to 1.0.

[0163] Figure 4 The inversion results were obtained using the method described in this application. The inversion results recover a low resistivity anomaly in the center of the survey area, whose horizontal location and burial depth range are consistent with... Figure 3 The actual anomaly shown is largely consistent with the real anomaly, indicating that the time-domain airborne electromagnetic 3D inversion based on POD-reduced order can identify the main low-resistivity structures in the real resistivity model. Due to the influence of the spatial resolution of the observation data, model regularization, and the non-uniqueness of 3D electromagnetic inversion, the boundary of the recovered anomaly is gentler than that of the real model and extends to some extent in the depth direction. Combined with the change in the normalized root mean square error (RMS), it can be seen that after switching from the full-order Gaussian-Newton iteration to the POD-reduced order iteration, the data fitting error continues to decrease, indicating that the full-order inversion parameter space and the POD-reduced order parameter space can achieve consistent transmission and maintain the ability to recover the target low-resistivity anomaly.

[0164] The above description is merely a preferred embodiment of this application and is not intended to limit this application. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of this application should be included within the protection scope of this application.

Claims

1. A time-domain airborne electromagnetic three-dimensional inversion method based on POD-reduced order, characterized in that, include: Determine the set of undersampled emission sources based on the spatial location of the emission sources in the test area; A local computational region is established for each undersampled emission source in the set of undersampled emission sources, and the three-dimensional electromagnetic response corresponding to each undersampled emission source is calculated. The three-dimensional electromagnetic response is converted into an undersampled time-domain predicted response for each measurement time channel, and the complete time-domain predicted response corresponding to each emission source position is iteratively reconstructed in the sparse transform domain. A full-order Gaussian-Newton inversion is performed. In each iteration, the full-order weighted residual and full-order weighted sensitivity matrix are calculated based on the predicted response in the complete time domain. The gradient of the objective function and the Gaussian-Newton approximation Hessian matrix in the full-order Gaussian-Newton inversion are constructed based on the full-order weighted residual and full-order weighted sensitivity matrix, and the full-order Gaussian-Newton equation is formed. The resistivity model update direction is obtained by solving the full-order Gaussian-Newton equation. The accepted full-order trial resistivity model is determined by full-order line search. When the preset order reduction switching condition is met, the last accepted full-order trial resistivity model before order reduction is selected as the POD source resistivity model. Multiple snapshots are constructed centered on the POD source resistivity model, forming a snapshot matrix. Singular value decomposition is then performed on the snapshot matrix to obtain the POD basis. The POD basis is mapped to the internal inversion parameter space to obtain a reduced-order basis; Using the reduced-order basis, the full-order weighted residual and the full-order weighted sensitivity matrix are projected onto the reduced-order space, and the reduced-order Gauss-Newton equation is solved to obtain the update direction of the reduced-order model; The order reduction model update direction is generated by searching the order reduction line to generate an order reduction trial resistivity model, and the complete time domain predicted response of the order reduction trial resistivity model is calculated; only the order reduction trial resistivity model that satisfies the objective function value descent condition is accepted as the current inversion model; When the updated current inversion model meets the preset termination condition, the final inversion model is output.

2. The time-domain airborne electromagnetic three-dimensional inversion method based on POD order reduction according to claim 1, characterized in that, The set of undersampled emission sources is determined based on their spatial location within the survey area, including: The Poisson disk sampling method is used to select the emission source positions in the two-dimensional spatial grid of the test area, such that the spatial distance between any two selected emission source positions is not less than a preset minimum spacing threshold, thus obtaining the undersampled emission source set.

3. The time-domain airborne electromagnetic three-dimensional inversion method based on POD order reduction according to claim 1, characterized in that, For each undersampled emission source in the set of undersampled emission sources, a local computational region is established, and the three-dimensional electromagnetic response corresponding to each undersampled emission source is calculated, including: Centered on the current undersampled emission source location, the local effective influence area is determined based on the skin depth; A local grid is established within the local effective influence area, and an extended cell layer is provided outward from the boundary of the local grid to absorb boundary reflections; Solving the secondary electric field diffusion equation yields the secondary electric field distribution at each node on the local grid. This distribution is then superimposed with the background field to obtain the three-dimensional electromagnetic response of the current undersampled emission source.

4. The time-domain airborne electromagnetic three-dimensional inversion method based on POD order reduction according to claim 1, characterized in that, The three-dimensional electromagnetic response is converted into an undersampled time-domain predicted response for each measurement time channel, and the complete time-domain predicted response corresponding to each emission source position is iteratively reconstructed in the sparse transform domain, including: The frequency domain response of each undersampled emission source is transformed into the frequency-time domain to obtain the undersampled time domain predicted response of each measurement time channel; For each measurement time channel, the undersampled time domain predicted response is arranged into a two-dimensional response matrix according to the spatial location of the emission source, the unsampled positions are set to zero, and a binary sampling position indication matrix of the same size as the two-dimensional response matrix is ​​established. Using curvelet transform as a sparse transform, the two-dimensional response matrix of each measurement time channel is iteratively reconstructed. Based on the binary sampling position indication matrix, the response values ​​of the sampled positions in the two-dimensional response matrix are forcibly replaced with the corresponding values ​​in the original undersampled time domain predicted response, thus completing the data reconstruction of each measurement time channel. After iterative convergence, the data reconstruction results are assembled in a predetermined order into complete time-domain predicted responses corresponding to each emission source location.

5. The time-domain airborne electromagnetic three-dimensional inversion method based on POD order reduction according to claim 1, characterized in that, In each iteration, the full-order weighted residual and full-order weighted sensitivity matrix are calculated based on the predicted response in the complete time domain. The gradient of the objective function and the Gauss-Newton approximation Hessian matrix in the full-order Gauss-Newton inversion are then constructed based on these matrices, forming the full-order Gauss-Newton equation. Solving the full-order Gauss-Newton equation yields the resistivity model update direction, including: The full-order weighted residuals are calculated based on the difference between the observed data and the predicted response in the complete time domain, combined with the data covariance matrix. Calculate the full-order weighted sensitivity of the predicted response in the full time domain to the resistivity model parameters under the current resistivity model; when the full-order weighted sensitivity is stored in the local window of the emission source, extract the corresponding row from the global network model space through index mapping, and combine them into a full-order weighted sensitivity matrix in the order of emission source, time step, measurement point, component and time channel; The gradient of the objective function is calculated from the full-order weighted residual and the full-order weighted sensitivity matrix. The Gauss-Newton approximation Hessian matrix is ​​constructed by multiplying the transpose of the full-order weighted sensitivity matrix with the full-order weighted sensitivity matrix, plus the product of the regularization parameter and the data covariance matrix. The linear equation system of the Gauss-Newton approximation Hessian matrix and the gradient of the objective function is solved to obtain the resistivity model update direction.

6. The time-domain airborne electromagnetic three-dimensional inversion method based on POD order reduction according to claim 1, characterized in that, Multiple snapshots are constructed centered on the POD source resistivity model, forming a snapshot matrix. Singular value decomposition is then performed on the snapshot matrix to obtain the POD basis, including: The spatial regions in the POD source resistivity model where the resistivity value deviates from the background resistivity by a preset value are marked as abnormal regions. Using the abnormal regions as perturbation objects, a set of snapshots is generated using multiple perturbation modes. The generated snapshots are arranged in columns to form a snapshot matrix. The snapshot matrix is ​​then processed to remove the mean and center the snapshots so that the mean of each snapshot is zero. Singular value decomposition is performed on the mean-removed snapshot matrix. The resulting singular values ​​are arranged in descending order. The square of each singular value represents the energy share of the corresponding mode in the snapshot matrix. The energy share of the first few singular values ​​is accumulated sequentially. When the proportion of the accumulated energy to the total energy first reaches the preset holding threshold, the number of modes contained in the current accumulation is used as the reduced dimension. Extract the left singular vectors corresponding to the selected reduced dimension from the singular value decomposition results to form the POD basis.

7. The time-domain airborne electromagnetic three-dimensional inversion method based on POD order reduction according to claim 1, characterized in that, Mapping the POD basis to the internal inversion parameter space yields a reduced-order basis, including: Based on the coordinate correspondence between regular grid cells and internal inversion cells, a constraint mapping from the regular grid index to the internal inversion cell index is established. The constraint mapping is used to filter and rearrange the POD base, remove redundant rows corresponding to non-inversion cells, and rearrange the rows corresponding to internal inversion cells according to the storage order of the internal inversion cells to obtain the POD base in the internal cell space. Using the square root of model covariance as the mapping operator, the mapping equation is solved for each column of the POD basis in the internal cell space to obtain the reduced-order basis of the internal inversion parameter space.

8. The time-domain airborne electromagnetic three-dimensional inversion method based on POD order reduction according to claim 7, characterized in that, Calculate the basis mapping consistency error, which is the relative deviation between the POD basis in the internal cell space and the reduced-order basis result of the model covariance square root operator acting on the internal inversion parameter space. When the basis mapping consistency error exceeds a preset threshold, stop entering the reduced-order iteration and revert to the full-order iteration.

9. The time-domain airborne electromagnetic three-dimensional inversion method based on POD order reduction according to claim 1, characterized in that, Using the reduced-order basis, the full-order weighted residuals and full-order weighted sensitivity matrices are projected into the reduced-order space, and the reduced-order Gauss-Newton equations are solved to obtain the update direction of the reduced-order model, including: Multiply the full-order weighted sensitivity matrix by the reduced-order basis to obtain the reduced-order weighted sensitivity matrix; Multiplying the reduced-order weighted sensitivity matrix by the full-order weighted residual yields the reduced-order negative gradient vector; The left-hand coefficient matrix of the reduced-order Gauss-Newton equation is formed by multiplying the transpose of the reduced-order weighted sensitivity matrix with the product of the regularization parameter and the transpose of the reduced-order basis with the product of the reduced-order basis; the reduced-order Gauss-Newton equation is constructed by using the reduced-order negative gradient vector as the right-hand term of the reduced-order equation. Solve the reduced-order Gauss-Newton equation to obtain the update direction of the reduced-order model.

10. The time-domain airborne electromagnetic three-dimensional inversion method based on POD order reduction according to claim 9, characterized in that, The reduced-order model update direction is generated by searching the reduced-order line to generate a reduced-order trial resistivity model, and the complete time-domain predicted response of the reduced-order trial resistivity model is calculated. Only reduced-order trial resistivity models that satisfy the objective function value descent condition are accepted as the current inversion model, including: With a step size of 1 as the initial value, a trial reduction coefficient is generated along the update direction of the reduction model. Based on the linear combination of the current internal inversion model, the reduction basis, and the trial reduction coefficient, a trial internal inversion model is generated. Based on the mapping relationship between the internal inversion parameters and the resistivity model parameters, the trial internal inversion model is converted into a reduced-order trial resistivity model. Calculate the complete time-domain predicted response of the reduced-order trial resistivity model; The objective function value is calculated based on the difference between the observed data and the predicted response in the complete time domain, and it is determined whether the objective function value satisfies the descent condition relative to the objective function value of the current inversion model. If the descent condition is met, the trial resistivity model is accepted as the current inversion model, and the corresponding full time domain prediction response, objective function value, and normalized root mean square error are updated. If the descent condition is not met, the trial order reduction coefficient, trial internal inversion model, and reduced-order trial resistivity model are regenerated after reducing the step size by a preset ratio. The full time domain prediction response is recalculated and re-evaluated until the descent condition is met or the maximum number of attempts is reached. If no reduced-order resistivity model satisfies the descent condition after reaching the maximum number of attempts, the current inversion model remains unchanged, and the search for the reduced-order line is terminated.

Images