Zhengyong Ren

3D direct current resistivity modeling with unstructured mesh by adaptive finite-element method

A new adaptive finite-element method for solving 3D direct-current resistivity modeling problems is presented. The method begins with an initial coarse mesh, which is then adaptively refined wherever a gradient-recovery-based a posteriori error estimator indicates that refinement is necessary. Then the problem is solved again on the new grid. The alternating solution and refinement steps continue until a given error criterion is satisfied. The method is demonstrated on two synthetic resistivity models with known analytical solutions, so the errors can be quantified. The applicability of the numerical method is illustrated on a 2D homogeneous model with a topographic valley. Numerical results show that this method is efficient and accurate for geometrically complex situations.

Fast 3‐D large‐scale gravity and magnetic modeling using unstructured grids and an adaptive multilevel fast multipole method

A novel fast and accurate algorithm is developed for large-scale 3D gravity and magnetic modeling problems. An unstructured grid discretization is used to approximate sources with arbitrary mass and magnetization distributions. A novel adaptive multi-level fast multipole (AMFM) method is developed to reduce the modeling time. An observation octree is constructed on a set of arbitrarily distributed observation sites, while a source octree is constructed on a source tetrahedral grid. A novel characteristic is the independence between the observation octree and the source octree, which simplifies the implementation of different survey configurations such as airborne and ground surveys. Two synthetic models, a cubic model and a half-space model with mountain-valley topography, are tested. As compared to analytical solutions of gravity and magnetic signals, excellent agreements of the solutions verify the accuracy of our AMFM algorithm. Finally, our AMFM method is used to calculate the terrain effect on an airborne gravity data set for a realistic topography model represented by a triangular surface retrieved from a digital elevation model. Using 16 threads, more than 5800 billion interactions between 1002001observation points and 5839830 tetrahedral elements are computed in 453.6 seconds. A traditional first-order Gaussian quadrature approach requires 3.77 days. Hence, our new AMFM algorithm not only can quickly compute the gravity and magnetic signals for complicated problems, but can also substantially accelerate the solution of 3D inversion problems.

A hybrid boundary element-finite element approach to modeling plane wave 3D electromagnetic induction responses in the Earth

A novel hybrid boundary element-finite element scheme which is accelerated by an adaptive multi-level fast multipole algorithm is presented to simulate 3D plane wave electromagnetic induction responses in the Earth. The remarkable advantages of this novel scheme are the complete removal of the volume discretization of the air space and the capability of simulating large-scale complicated geo-electromagnetic induction problems. To achieve this goal, first the Galerkin edge-based finite-element method (FEM) using unstructured meshes is adopted to solve the electric field differential equation in the heterogeneous Earth, where arbitrary distributions of conductivity, magnetic permeability and dielectric permittivity are allowed for. Second, the point collocation boundary-element method (BEM) is used to solve a surface integral formula in terms of the reduced electrical vector potential on the arbitrarily shaped air-Earth interface. Third, to avoid explicit storage of the system matrix arising from large-scale problems and to reduce the horrendous time complexity of the product of the system matrix with an initial vector of unknowns, the adaptive multilevel fast multipole method is applied. This leads to a matrix-free form suitable for the application of iterative solvers. Furthermore, a highly sparse problem-dependent preconditioner is developed to significantly reduce the number of iterations used by the iterative solvers.The efficacy of the presented hybrid scheme is verified on two synthetic examples against different numerical techniques such as goal-oriented adaptive finite-element methods. Numerical experiments show that at low frequencies, where the quasi-static approximation is applicable, standard FEM methods prove to be superior to our hybrid BEM-FEM solutions in terms of computational time, because the FEM method requires only a coarse discretization of the air domain and offers an advantageous sparsity of the system matrix. At radio-magnetotelluric frequencies of a few hundred kHz, the hybrid BEM-FEM scheme outperforms the FEM method, because it avoids explicit storage of the system matrices as well as dense volume discretization of the air domain required by FEM methods at high frequencies. In summary, to the best of our knowledge, this study is the first attempt at completely removing the air space for large scale complicated electromagnetic induction modeling in the Earth.

A new extrapolation cascadic multigrid method for three dimensional elliptic boundary value problems

In this paper, we develop a new extrapolation cascadic multi grid method (ECMGjcg), which makes it possible to solve 3D elliptic boundary value problems on r ectangular domains of over 100 million unknowns on a desktop computer in minutes. First, by combining Richar dson extrapolation and tri-quadratic Serendipity interpolation techniques, we introduce a new extrapolatio n formula to provide a good initial guess for the iterative solution on the next finer grid, which is a third order approxi mation to the finite element (FE) solution. And the resulting large sparse linear system from the FE discretiza tion is then solved by the Jacobi-preconditioned Conjugate Gradient (JCG) method. Additionally, instead of performin g a fixed number of iterations as used in the most of cascadic multigrid method (CMG) literature, a relative res idual stopping criterion is used in our iterative solvers, which enables us to obtain conveniently the numerical solut ion with the desired accuracy. Moreover, a simple Richardson extrapolation is used to cheaply get a fourth ord er accurate solution on the entire fine grid from two second order accurate solutions on two di fferent scale grids. Test results from three di fferent problems with smooth and singular solutions are reported to show that ECMG jcg has much better e fficiency compared to the classical Vcycle and W-cycle multigrid methods. Since the initial gues s for the iterative solution is a quite good approximation to the FE solution, numerical results show that only few numb er of iterations are required on the finest grid for ECMGjcg with an appropriate tolerance of the relative residual to ac hieve full second order accuracy, which is particularly important when solving large systems of equat ions and can greatly reduce the computational cost. It should be pointed out that when the tolerance becomes more cr uel, ECMGjcg still needs only few iterations to obtain fourth order extrapolated solution on each grid, except on t he finest grid. Finally, we present the reason why our ECMG algorithms are so highly e fficient for solving these elliptic problems.In this paper, we develop a new extrapolation cascadic multigrid method, which makes it possible to solve three dimensional elliptic boundary value problems with over 100 million unknowns on a desktop computer in half a minute. First, by combining Richardson extrapolation and quadratic finite element (FE) interpolation for the numerical solutions on two-level of grids (current and previous grids), we provide a quite good initial guess for the iterative solution on the next finer grid, which is a third-order approximation to the FE solution. And the resulting large linear system from the FE discretization is then solved by the Jacobi-preconditioned conjugate gradient (JCG) method with the obtained initial guess. Additionally, instead of performing a fixed number of iterations as used in existing cascadic multigrid methods, a relative residual tolerance is introduced in the JCG solver, which enables us to obtain conveniently the numerical solution with the desired accuracy. Moreover, a simple method based on the midpoint extrapolation formula is proposed to achieve higher-order accuracy on the finest grid cheaply and directly. Test results from four examples including two smooth problems with both constant and variable coefficients, an H3-regular problem as well as an anisotropic problem are reported to show that the proposed method has much better efficiency compared to the classical V-cycle and W-cycle multigrid methods. Finally, we present the reason why our method is highly efficient for solving these elliptic problems.

Gravity Gradient Tensor of Arbitrary 3D Polyhedral Bodies with up to Third-Order Polynomial Horizontal and Vertical Mass Contrasts

During the last 20 years, geophysicists have developed great interest in using gravity gradient tensor signals to study bodies of anomalous density in the Earth. Deriving exact solutions of the gravity gradient tensor signals has become a dominating task in exploration geophysics or geodetic fields. In this study, we developed a compact and simple framework to derive exact solutions of gravity gradient tensor measurements for polyhedral bodies, in which the density contrast is represented by a general polynomial function. The polynomial mass contrast can continuously vary in both horizontal and vertical directions. In our framework, the original three-dimensional volume integral of gravity gradient tensor signals is transformed into a set of one-dimensional line integrals along edges of the polyhedral body by sequentially invoking the volume and surface gradient (divergence) theorems. In terms of an orthogonal local coordinate system defined on these edges, exact solutions are derived for these line integrals. We successfully derived a set of unified exact solutions of gravity gradient tensors for constant, linear, quadratic and cubic polynomial orders. The exact solutions for constant and linear cases cover all previously published vertex-type exact solutions of the gravity gradient tensor for a polygonal body, though the associated algorithms may differ in numerical stability. In addition, to our best knowledge, it is the first time that exact solutions of gravity gradient tensor signals are derived for a polyhedral body with a polynomial mass contrast of order higher than one (that is quadratic and cubic orders). Three synthetic models (a prismatic body with depth-dependent density contrasts, an irregular polyhedron with linear density contrast and a tetrahedral body with horizontally and vertically varying density contrasts) are used to verify the correctness and the efficiency of our newly developed closed-form solutions. Excellent agreements are obtained between our solutions and other published exact solutions. In addition, stability tests are performed to demonstrate that our exact solutions can safely be used to detect shallow subsurface targets.

A finite-volume approach for 2D magnetotellurics modeling with arbitrary topographies

A novel finite-volume approach for complicated 2D magnetotellurics (MT) problems with arbitrarily surface topography is presented. An edge-surface integral balance equation is derived by employing a conservation law on the generalized 2D MT boundary value problem. A triangular grid is used to discretize the 2D conductivity model so that we can deal with arbitrarily complex cases with surface topography. The node-centered finite-volume algorithm is used to derive the final system of linear equations on a dual mesh of the triangular grid, which is solved by a robust direct solver. Three synthetic models verify the accuracy of the presented finite-volume algorithm and its capability of dealing with surface topography.

Analytical Formulas for Underwater and Aerial Object Localization by Gravitational Field and Gravitational Gradient Tensor

Object localization techniques have significant applications in civil fields and safety problems. A novel analytical formula is developed for accurate underwater and aerial object real-time localization by combining gravitational field and horizontal gravitational gradient anomalies. The proposed method enhances the accuracy of object localization and its excess mass estimation; it also effectively avoids the possible numerical instability and the singularity in the previous works. Finally, a synthetic underwater object navigation model was adopted to verify its performance. The results show that our newly developed method is more practical than existing methods.

Space-time array difference magnetotelluric method

Summary form only given. A space-time array difference magnetotelluric method is proposed in order to suppress the correlated noise and obtain reliable plane wave impedance data in strong interference region. First, we established a model of electromagnetic prospecting system with multi-input and multi-output. Both natural and man-made electromagnetic field is considered to incidence to earth simultaneously, constitute a space-time array input. Many stations are carried out to simultaneously observe the total responses of all the inputs on the earths surface, constitute a space-time array output. A set of quadratic space-time equations is obtained by resolving the input-output relationship of the linear time invariant system. Second, we presented a four step scheme to solve the space-time equations within complex input environment. Step 1, three types of space-time array data matrixes are obtained by reasonable observation design in the field, including the target data matrix constructed by data from all the stations in the survey area, the natural field data matrixes constructed by data from the remote reference stations, and the man-made electromagnetic field data matrix constructed by data from the horizontal magnetic field differential signals of the survey stations. Step 2, by using the principal component analysis method, the polarization parameters of the natural sources and man-made noise sources are extracted from the natural field data matrix and man-made electromagnetic field data matrix respectively. Step 3, the system responses of all the observation channels corresponding to each sources are estimated from the target data matrix using linear regression method. Step4, impedance tensors, apparent resistivities and phases corresponding to different sources are estimated from the system responses, so as to achieve the goal of signal noise separation based on the input end. Last, numerical simulation and practical experiment are carried out to evaluate the proposed method. The results indicate that our method is better than the conventional magnetotelluric method for processing the data contaminated with correlated noise, and can effectively separate the responses of natural and man-made electromagnetic field. A larger size of the space-time array can obtain more reasonable results when the noise environment is complex in the survey area.

2.5-D DC resistivity forward modeling and inversion by finite element-infinite element coupled method

To reduce the numerical errors arising from the improper enforcement of the artificial boundary conditions on the distant surface that encloses the underground part of the subsurface, we present a finite-element–infinite-element coupled method to significantly reduce the computation time and memory cost in the 2.5D direct-current resistivity inversion. We first present the boundary value problem of the secondary potential. Then, a new type of infinite element is analysed and applied to replace the conventionally used mixed boundary condition on the distant boundary. In the internal domain, a standard finite-element method is used to derive the final system of linear equations. With a novel shape function for infinite elements at the subsurface boundary, the final system matrix is sparse, symmetric, and independent of source electrodes. Through lower upper decomposition, the multi-pole potentials can be swiftly obtained by simple back-substitutions. We embed the newly developed forward solution to the inversion procedure. To compute the sensitivity matrix, we adopt the efficient adjoint equation approach to further reduce the computation cost. Finally, several synthetic examples are tested to show the efficiency of inversion.

2.5D direct-current resistivity forward modelling and inversion by finite-element–infinite-element coupled method

To reduce the numerical errors arising from the improper enforcement of the artificial boundary conditions on the distant surface that encloses the underground part of the subsurface, we present a finite-element–infinite-element coupled method to significantly reduce the computation time and memory cost in the 2.5D direct-current resistivity inversion. We first present the boundary value problem of the secondary potential. Then, a new type of infinite element is analysed and applied to replace the conventionally used mixed boundary condition on the distant boundary. In the internal domain, a standard finite-element method is used to derive the final system of linear equations. With a novel shape function for infinite elements at the subsurface boundary, the final system matrix is sparse, symmetric, and independent of source electrodes. Through lower upper decomposition, the multi-pole potentials can be swiftly obtained by simple back-substitutions. We embed the newly developed forward solution to the inversion procedure. To compute the sensitivity matrix, we adopt the efficient adjoint equation approach to further reduce the computation cost. Finally, several synthetic examples are tested to show the efficiency of inversion.