任政勇

This page presents the research progress in 2022 from Zhengyong's group.  

1.  Geo-electromagentic induction modeling and inversion

1.1 复杂地质模型建模

当前，基于非结构化四面体网格的正反演技术具有逼近任意起伏地形和特定地质结构的能力。但任意复杂地质模型的四面体网格建立非常困难，制约了该技术在实际勘探中的应用。针对这一难题，本文提出了一套复杂地质-地球物理模型的建模流程，研发了一系列建模代码，具有融合数字高程数据、地质剖面图和钻孔数据等地质地理信息的能力，能够高效建立面向有限元模拟计算的高质量复杂地球物理模型。三维地质模型是表达各地质要素之间接触关系的可视化数字模型，是地质研究的基础。而三维地球物理模型是以网格单元为基础表示地下几何界面的一种定量化计算模型，其本质是网格节点、拓扑关系以及物性值的规则数值序列。尽管二者存在区别，但是其几何界面都是由二维的离散网格构建的。因此，可以利用三维网格剖分技术将三维地质模型转换为三维地球物理模型。本文提出了一套复杂地质-地球物理模型建模流程，该流程从多源数据融合建模方法出发，以SKUA-Gocad建模平台和Facetmodeller软件来建立三维地质模型，紧接着基于约束2D/3D-Delaunay网格剖分算法，对三维地质模型进行界面三角网格质量优化、局部网格加密以及相交界面网格约束重采样操作，从而建立一体化的三维体积模型，最后利用Tetgen网格剖分器生成用于有限元计算的高质量四面体网格。此流程建模效率较高，适用性较广且模型的修改优化和质量控制较为容易。为了验证本文所提出的建模流程的适用性，我们以铜陵冬瓜山铜矿区为例，建立了高质量的复杂地球物理模型，并利用之前推导的基于四面体的解析解公式和有限元算法计算了该模型的重磁、大地电磁（MT）以及可控源电磁（CSEM）响应。一系列的正演、反演计算结果表明本文提出的建模流程具有较高的适用性，建立的地球物理模型网格质量较好，能够满足不同算法的稳定计算，有广泛的潜在应用价值。本文提出的建模新技术，不仅可以用于消去矿山已知浅部地质体的干扰，而且可为深部数据反演提供高质量的参考模型，从而提高深部勘探的效果。

（任政勇，王祥，铜陵冬瓜山铜矿区复杂地球物理模型构建及其应用，地球物理学报，2022）

1.2  起伏地形情况下三维大地电磁自适应反演

当前，三维大地电磁反演常采用固定反演网格。过稀的反演网格无法逼近复杂的地下电导率分布,还可能导致反演结果无法收敛而过密的反演网格不仅将加大反演的计算成本，还会增加反演的非唯一性。针对上述问题，本文提出了一种基于灵敏度矩阵和模型梯度共同驱动的三维大地电磁自适应反演策略，以实现反演网格的自适应调整。首先，基于非结构化有限元技术实现了对任意复杂带地形三维地电模型的精确正演；其次，生成嵌套的正、反演四面体网格，实现了正反演网格的解耦；建立基于光滑约束的Thikhonov正则化目标函数，采用L-BFGS最优化算法、基于伴随原理的目标函数梯度求解策略以及正则化因子“冷却”技术实现了对目标函数的最优化求解。然后，提出了一种基于灵敏度矩阵和模型梯度双重约束的反演网格调整策略，实现了最优化过程中的反演网格自适应调整。最后，通过对多个理论模型的合成数据和实测数据的反演验证了该自适应反演算法的可靠性、稳定性和实用性。

(Huang Chen, Zhengyong Ren*, Jingtian Tang, et.al. JGR-solid earth, 2022, under preparing.)

1.3 Marine CSEM modeling

We present a goal-oriented adaptive finite-element algorithm for accurately modelling marine controlled-source electromagnetic responses in 3-D media with general electrical anisotropy. We formulate the primal boundary value problem in terms of the total electric field for general applications. Following the goal-oriented adaptivity concept, the dual problem is derived from a functional designed to measure the data errors of interest in light of the data quality. We approximate the solutions to the primal and dual formulations using the edge finite-element method on tetrahedral grids for a flexible treatment of complex geological settings and survey geometries. To control the mesh adaptation, we develop a reliable residual-type a posteriori error estimation that takes account of the volumetric residual and the numerical discontinuity of the normal component of the electrical current density and that of the tangential component of the magnetic field with respect to non-smooth and anisotropic coefficients.We demonstrate the proposed modelling solver on 1-D MCSEM scenarios with varying degrees of electrical anisotropy. The comparison with goal-oriented adaptivity results obtained from other three commonly used error indicators shows that our approach is robust in dealing with both moderate and strong electrical anisotropy. After that, we constructed a 3-D hydrocarbon-bearing reservoir model with slope seafloor topography and tilted transverse isotropy in the background to examine our algorithm for the case of multiple sources. Finally, we implement a sensitivity analysis procedure to evaluate the resolution of the electrical anisotropy. The quantitative results indicate limitations and preferences of conventional MCSEM data in resolving anisotropic models, providing fundamental insights for inversion based data interpretation. 

(Feiyan Wang, Zhengyong Ren*, Lihong Zhao, A goal-oriented adaptive finite element approach for 3D marine controlled-source electromagnetic problems with general electrical anisotropy, Geophysical Journal International, https://doi.org/10.1093/gji/ggab485, 2022)

1.4 Land CSEM modeling

The CSEM method, as an essential method in exploration geophysics, which has been widely used for different geological and geophysical problems, such as mineral deposits exploration, offshore hydrocarbon exploration and reservoir monitoring. To efficiently interpret or invert controlled-source electromagnetic (CSEM) field data which are more and more recorded at surveying areas with complex geological environment and arbitrary topography, three-dimensional (3D) CSEM forward modeling softwares with capabilities of fast solving large-scale problems, offering accurate electromagnetic responses for complex geo-electrical models and particularly easily being incorporated into inversion are required. In this study, we developed a parallel goal-oriented adaptive mesh refinement (AMR) finite-element approach for frequency-domain 3D CSEM forward modelling with hierarchical tetrahedral grids, which can offer accurate electromagnetic responses for large-scale complex models and efficiently serve for inversion.

(Zhengguang Liu, Zhengyong Ren*,  Hongbo Yao, Jingtian Tang, etal. A parallel adaptive finite-element approach for three-dimensional realistic controlled-source electromagnetic problems using hierarchical tetrahedral grids. Geophysics, 2022, to be submitted.)

1.5 Land CSEM inversion

Controlled-source electromagnetic (CSEM) method is crucial for detecting and locating underground anomalies and structures. However, it is challenging to interpret the field data with multi-frequency via 3D CSEM inversion. To fully excavate and utilize the valuable information of CSEM data at different frequencies, we propose an efficient algorithm for 3D multifrequency CSEM (MFCSEM) inversion based on rational Krylov (RK) subspace. Within the framework of our algorithm, we first use the three-term Lanczos recursion to construct the RK basis matrix quickly; thus, the fast MFCSEM forward modeling can be realized via the RK approximation. Subsequently, we present a novel cyclic projection and correction (CPC) algorithm to solve the MFCSEM adjoint forward problems. Finally, the nonlinear conjugate gradient (NLCG) method is adopted to seek solutions to the nonlinear inverse problem. We demonstrate the excellent performance of our algorithm by synthetic and field data sets. The inversion results show that our algorithm is computationally efficient, resulting in considerable speedup compared with the conventional method. Our algorithm provides a new idea that would significantly improve the efficiency of MFCSEM inversion.

(Liu J R, Ren, ZY, etc, IEEE Transactions on Geoscience and Remote Sensing, 2022, under revision)

1.6  Global EM induction modeling

We present a multi-resolution finite-element approach for three-dimensional (3D) electromagnetic (EM) induction modeling in spherical Earth excited by external current sources. Firstly, the secondary electric field is utilized to formulate the boundary value problem so that both magnetospheric and ionospheric current sources are naturally considered. Secondly, the conforming tetrahedral grids with multi-resolutions or multiple sizes are used to efficiently approximate the heterogeneous crust and mantle, especially for ocean-continent interfaces, so that the local ocean effects at coastal and island observatories can be accurately simulated. Furthermore, a parallel goal-oriented hp-adaptive finite-element method with Nedelec vector elements is employed to guarantee the accuracy of final solutions for arbitrary 3D conductivity distributions. Finally, two synthetic models are used to verify the accuracy and efficiency of our newly developed forward modeling solver. Results show that accurate solutions can be obtained for problems with several millions to hundred millions unknowns in a few minutes on a cluster using 128 cores, which is comparable to the efficiency of the mainstream integral equation forward solver. We apply this approach to correct the near-surface ocean effects for several unused Chinese coastal observatories by performing multi-resolution 3D modeling. The corrected geomagnetic data are inverted for the subsurface layered mantle conductivity structures. The conductivity model beneath southeast China is more resistive than that beneath northeast China by more than half an order of magnitude. By comparing the inverse models with the latest laboratory conductivity-depth profiles, the estimated transition zone water content is less than 0.01wt% beneath southeast China irrespective of which laboratory data is used. Considering the low velocity anomalies in this region, which suggest high temperature structures, less water is expected to explain the low conductivity anomalies. We therefore infer that the mantle transition zone beneath southeast China is dry.

( Yao H,  Ren Z*,  A multi-resolution finite-element approach for global electromagnetic induction modeling with application to Chinese coastal geomagnetic observatory studies, JGR-SOLID EARTH,  under preparing.)

1.7 Earth mantle conductivity inversion

China has a long coastline which is surrounded by active subduction zones. The subducted Pacific and Philippine Sea plates have played important roles in the tectonic evolution of eastern China. The water in the subducting oceanic slabs can be carried into the mantle transition zone and the lower mantle. A small amount of water can significantly affect the physical and chemical properties of mantle minerals, such as the plastic deformation and the melting temperature. Therefore, constraining the distribution and amount of water in the mantle beneath China plays an important role in understanding the dynamics and evolution of the Earth. Chinese coastal geomagnetic observatory data provides valuable data source to study the wa ter content of subduction zone by using global EM induction method. Recently, the coastal geomagnetic data in northeast China have been used to constraint the electrical conductivity and water content. However, several coastal geomagnetic observatories in southeast China (Quanzhou and Yongning) are not used due to the lack of a multi-resolution forward modeling solver to accurately account for the local ocean effects (Kuvshinov et al. 2002). Having validated the accuracy and efficiency of our solver on synthetic models, we now apply our multi-resolution finite element approach to guide the interpretation of Chinese coastal observatory data. This application aims to answer the following questions: (1) How the ocean affect the responses from coastal geomagnetic observatories in China? (2) How dense the resolution of surafce layer is needed to accurately model the local ocean effects? (3) Will inaccurate correction for the ocean effects lead to large difference among inverted conductivity profiles? (4) Is the transition zone wet or dry beneath southeast China?

( Yao H,  Ren Z*,  A multi-resolution finite-element approach for global electromagnetic induction modeling with application to Chinese coastal geomagnetic observatory studies, JGR-SOLID EARTH,  under preparing.)

Mantle conductivity imaging is one of scientific goals of the forthcoming Macau's scientific satellite. To achieve this goal, we develop a data analysis and inversion scheme for satellite magnetic data to probe global one-dimensional (1D) mantle conductivity structures. Using this scheme, we present a new global mantle conductivity model derived from over 8 years of Swarm satellite magnetic data. This scheme can be used to process the forthcoming Macau's scientific satellite magnetic data.

(Yao H, Ren Z*,  Pan K, Tang J, Zhang K. A global mantle conductivity model derived from 8 years of Swarm satellite magnetic data. Earth and Planetary Physics, 2022, under review.)

1.8 磁暴产生的瞬时场模拟

 地球电磁感应法能获得地球深部电导率的分布，从而被广泛应用于探索地球内部结构、物质组分和热力学状态等重要信息，对正确认识地球演化过程及动力学机制等具有十分重要的意义。传统的地球电磁感应研究一般在频率域进行，然而，观测数据基本上为时间域电磁场信号，因此，直接分析时间域地球电磁感应数据具有先天优势。本文开发了一套基于时间域有限元法的地球电磁感应三维正演求解器，实现了高精度地计算地球外部磁层环电流产生的感应电磁场时间序列。首先，利用磁场矢量势建立了时间域地球电磁感应控制方程，结合磁层环形电流源的物理属性建立了边界条件和初始条件；然后，利用四面体矢量有限元和隐式后退欧拉公式，分别实现磁场矢量势的空间域和时间域离散，获得了不同时刻的实系数大型线性方程组，借助于MUMPS高性能并行直接求解器，快速高精度求解不同时刻的磁场矢量势和感应磁场。最后，利用理论模型分别在频率域和时间域中验证了本文算法的正确性；并利用Dst观测指数建立的真实激励场源，研究了中国和澳大利亚下方地幔转换带深度410km-670km范围内高导体的敏感性和可探测性，表明这些高导体在时间域内产生了明显的异常，证明了在时间域分析地球电磁感应数据的优越性能。

（任政勇，杨聪等，中国科学，2022，to be submitted）

2. Gravity and magnetic data modeling and inversion

2.1 Magnetic closed-form solution 

Magnetization is a natural property of magnetic materials which is widely used to study the ancient climate change whose footprint was recorded in rocks by palaeomagnetic data, uncover the location of unexploded ordnance, mineral deposits and gas and oil resources, and image regional crust and planet structure. An essential but very challenging task is to find analytical solutions of induced magnetic signals of general magnetic targets. In this study, for the first time, we demonstrate the existence of such solutions and explicitly give the analytical expressions of magnetic potential, magnetic field and magnetic gradient tensor for magnetic materials with polyhedral shapes and polynomial magnetization vectors. The order of magnetization polynomials can vary from linear order to arbitrary integer orders along both horizontal and vertical directions. Several synthetic and realistic iron deposit models validate the accuracy and eciency of newly developed analytical solutions. Due to its inherent importance, our findings would have an immediate impact on magnetic explorations in earth science and have potentials in dierent applications.

(Zhengyong Ren, Yuanlei Zhang, Yiyuan Zhong, Kejia Pan, QihongWu and Jingtian Tan, Analytical anomalies of polyhedral bodies with polynomial magnetization, NCS, submitted, 2022) 

2.2 重力场自适应有限元模拟

地球重力场作为地下介质密度场源产生的天然物理场，在地球科学中有丰富的研究历史和重要价值。重力场数据包含了地下密度结构的特征响应，被广泛用于矿产资源勘查、工程地质等领域。随着地面、航空和卫星重力技术的发展，海量重力场数据反演与解释是如今研究的热点，而制约反演发展的重要因素是正演的效率与精度。正演主要分为两个步骤：网格建模和计算，这两个部分共同决定了正演的效果，本文从这两方面入手分别进行了研究。为了对真实大尺度模型进行网格建模与剖分，本文基于真实地球圈层模型分别构建了球体三棱柱-四面体网格、球体四面体网格和球体六面体网格，并分析了它们的适用性和优缺点，同时还考虑了地球圈层模型中的局部模型，提出了构建大规模区域网格的思路和方案，并以青藏高原和台湾岛为例进行了构建。在此基础上，为了让正演达到效率与精度的平衡，本文基于非结构化网格提出了重力场及梯度的h-型自适应有限元三维正演方法。首先从重力场满足的泊松方程出发，推导了目前重力场已有的有限元方法和不同的边界条件，选择了最佳的方法基于伽辽金加权残差法将边值问题转换为求解大型稀疏线性方程组，从而得到有限元数值解。为了解决网格分布不合理的情况，引入了h-型自适应加密算法进行网格加密，通过对理论立方体模型和真实小行星Bennu模型的测试，验证了本文网格建模和算法的精度与效率。 

（任政勇，岳国旋，基于四面体自适应有限元的重力场及梯度计算方法, Geophysics, 2022, to be submitted.）

2.3 Gravity inversion using satelite data. 

We develop a novel gravity inversion algorithm in spherical coordinates based on adaptive inversion mesh refinement and multi-physical parameter constraints. The inversion mesh is discretized into tesseroids (spherical prisms) to take the curvature of the Earth into account. To reduce the number of unknowns and computational cost, the inversion mesh is adaptively refined according to the spatial variation of parameters at each iteration. Wavelet compression is used to further reduce the computational requirement. Besides, to alleviate the ambiguity of inversion, the algorithm is capable of incorporating a priori models from other geophysical methods via the cross-gradient coupling and the direct parameter coupling. For the cross-gradient coupling, the vector product of spatial gradients of two model parameters is included in the objective function. For the direct parameter coupling, the a priori model is converted to a density model using the relation between two physical properties, which is then used as the initial and reference model in the inversion. A synthetic example is presented to demonstrate the effectiveness of our inversion method. Finally, we invert the gravity data over the Tibetan plateau to obtain the density variations in the upper mantle, with a global S-wave tomographic model used as a constraint. The inversion model shows that the mantle lithosphere beneath the Himalayan collision zone varies from west to east. Low-density anomalies are observed beneath the southern and the northern Tibetan plateau, which are consistent with published low-velocity anomalies and magmatism distributions, possibly indicating two asthenospheric upwellings.

(Zhong Y,  Ren Z*, Tang J.T,  Constrained gravity inversion with adaptive inversion grid refinement in spherical coordinates and its application to mantle structure beneath Tibetan plateau,  JGR-SOLID EARTH, 2022, https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2021JB022916)

2.4 Gravity inversion using deposit data

In large-scale inversion problem, discretization of the inversion region may result in a tremendous number of elements. In addition, gravity inversion suffers from severe ambiguity and the increased number of unknowns might further aggravate the non-uniqueness. We use an adaptively refined mesh in gravity inversion to reduce the number of inversion unknowns. A priori information can be exploited to reduce the non-uniqueness of the inversion.

(Zhong Y,  Ren Z*, Tang J.T,  Adaptive gravity inversion using satalitte gravity data sets,  Geophysics, 2022, to be submitted.)