时间域航空电磁各向异性大地三维自适应有限元正演研究

各向异性介质模型电性结构复杂,如何进行合理的网格剖分成为获得高精度正演结果的关键,为此本文开展时间域航空电磁各向异性大地三维自适应有限元正演算法研究.通过结合非结构时间域有限元算法和自适应网格优化技术,实现各向异性介质条件下三维时间域航空电磁自适应正演.考虑到时间域航空电磁响应随时间的衰减特性,为了综合评价不同时刻的后验误差,本文将时间作为加权因子,调整各个时刻后验误差的相对权重,进而实现对浅部和深部网格的同步优化.通过与一维解析结果进行对比验证了本文算法的可靠性.数值实验结果显示电导率各向异性对自适应网格影响严重,其最大主轴电导率的数值及其分布特征直接决定了网格加密效果.此外,各向异性对时间域航空电磁三分量响应的分布形态和异常幅值也会产生严重影响,利用全域视电阻率极性图,可以很好地识别各向异性主轴方向.

     

从瞬变电磁扩散场到拟地震波场的全时域反变换算法

将瞬变电磁满足的扩散方程转变为波动方程,然后利用地震类成像方法实现瞬变电磁虚拟波场成像,是实现瞬变电磁三维反演的有效手段之一.为了实现由扩散场到虚拟波场的转换,文中采用预条件正则化共轭梯度法求解波场反变换问题.首先,对几种离散方式进行比较,采用条件数最小的离散方式进行离散;然后选择最优的正则化参数,并利用超松弛预条件技术对系数矩阵进行预条件处理;最后,利用共轭梯度法进行迭代求解.超松弛预条件有效降低了系数矩阵的条件数,正则化方法使得反变换得到的波场稳定、可靠,共轭梯度法能够保证计算快速收敛.将反变换结果与已知虚拟波场函数对比,证明算法稳定、可信.将文中算法结果与前人研究结果进行对比,说明方法效果.通过实测数据的波场变换处理给出了文中方法的实际应用效果.结合反变换算法,对不同参数模型进行分析,总结了虚拟波场在色散介质中的传播规律.     

均匀半空间表面大定源瞬变电磁响应的快速算法

大定源是地面瞬变电磁观测的主要方式之一,该方式在大深度、高密度的面积测量时具有明显的优势,但当接收点偏离发射框中心时,由于场源的非对称性(即框边影响),数据处理和解释比较困难,特别是全程瞬变响应的精确计算相当耗时,本文介绍一种数值算法,它既能实现快速计算,又能满足精度要求,该算法通过对瞬变场垂直分量(bz)及其时间变化率((a)bz/(a)t)的核函数Y(Z)和Y＇(Z)表现特性的研究,以参数Z把整个瞬变过程分为早期阶段(Z→0),中期阶段和晚期阶段(Z→∞).计算全程响应时,早期和晚期阶段分别采用Y(Z)和Y＇(Z)的渐近表达式;对中期阶段,内层积分(即误差函数erf)采用有理Chebyshev渐近展开式,外层积分采用Romberg数值积分法.理论模型计算表明,利用该算法可以快速计算空间任意点(除发射边框)的全程响应核函数Y(Z)和Y'(Z).当测点到边框的距离大于边长的25％时,计算速度比常规数值积分算法快7倍;其它测点处计算速度比常规数值积分算法快4倍.全程时段的相对误差＜0.0002％.     

A wetting and drying algorithm with a combined pressure/free-surface formulation for non-hydrostatic models

A wetting and drying method for free-surface problems for the three-dimensional, non-hydrostatic Navier–Stokes equations is proposed. The key idea is to use a horizontally fixed mesh and to apply different boundary conditions on the free-surface in wet and dry zones. In wet areas a combined pressure/free-surface kinematic boundary condition is applied, while in dry areas a positive water level and a no-normal flow boundary condition are enforced. In addition, vertical mesh movement is performed to accurately represent the free-surface motion. Non-physical flow in the remaining thin layer in dry areas is naturally prevented if a Manning–Strickler bottom drag is used. The treatment of the wetting and drying processes applied through the boundary condition yields great flexibility to the discretisation used. Specifically, a fully unstructured mesh with any finite element choice and implicit time discretisation method can be applied. The resulting method is mass conservative, stable and accurate. It is implemented within Fluidity-ICOM [1] and verified against several idealized test cases and a laboratory experiment of the Okushiri tsunami.     

Local transmitting boundaries for transient elastic analysis

The aim of this paper is to investigate and develop alternative methods of analyzing problems in dynamic soil–structure-interaction (SSI). The interaction means that the amplitude of structural response is effected by additional energy dissipation through radiation and material damping in the soil. The surrounding soft soil behaves as a natural damper for a massive and stiff structure supported or embedded in it. The main focus is the major difficulty posed by such an analysis — the phenomena of waves that radiate outward from the excited structures towards infinity. In numerical calculations only a finite region of the foundation medium is analyzed and something is done to prevent the outgoing radiation waves from reflecting at the boundary region.Development of a simple and efficient finite element (FE) procedure for the solution directly in the time domain of transient SSI problems is the main concern. The central feature of the procedure is local absorbing boundaries used to render the computational domain finite. These boundaries are local in both time and space and are completely defined by a pair of symmetric stiffness and damping matrices. As the effort for implementing them is the same as for the impedance boundary condition (BC) considering the angle of incidence, standard assembly procedure can be used. Due to the local nature they also preserve the overall structure of the global equations of motion. Even though the focus is in the time domain the same equations of motions can be used to determine the solution under time-harmonic excitation directly in the frequency domain. Explicit formulae for the element matrices are included in the paper and numerical examples for transient radiation model problems to illustrate the validity and accuracy of the new procedures, are given.     

A least-squares penalty method algorithm for inverse problems of steady-state aquifer models

Based on the generalized Gauss–Newton method, a new algorithm to minimize the objective function of the penalty method in (Bentley LR. Adv Wat Res 1993;14:137–48) for inverse problems of steady-state aquifer models is proposed. Through detailed analysis of the “built-in” but irregular weighting effects of the coefficient matrix on the residuals on the discrete governing equations, a so-called scaling matrix is introduced to improve the great irregular weighting effects of these residuals adaptively in every Gauss–Newton iteration. Numerical results demonstrate that if the scaling matrix equals the identity matrix (i.e., the irregular weighting effects of the coefficient matrix are not balanced), our algorithm does not perform well, e.g., the computation cost is higher than that of the traditional method, and what is worse is the calculations fail to converge for some initial values of the unknown parameters. This poor situation takes a favourable turn dramatically if the scaling matrix is slightly improved and a simple preconditioning technique is adopted: For naturally chosen simple diagonal forms of the scaling matrix and the preconditioner, the method performs well and gives accurate results with low computational cost just like the traditional methods, and improvements are obtained on: (1) widening the range of the initial values of the unknown parameters within which the minimizing iterations can converge, (2) reducing the computational cost in every Gauss–Newton iteration, (3) improving the irregular weighting effects of the coefficient matrix of the discrete governing equations. Consequently, the example inverse problem in Bentley (loc. cit.) is solved with the same accuracy, less computational effort and without the regularization term containing prior information on the unknown parameters. Moreover, numerical example shows that this method can solve the inverse problem of the quasilinear Boussinesq equation almost as fast as the linear one.In every Gauss–Newton iteration of our algorithm, one needs to solve a linear least-squares system about the corrections of both the parameters and the groundwater heads on all the discrete nodes only once. In comparison, every Gauss–Newton iteration of the traditional method has to solve the discrete governing equations as many times as one plus the number of unknown parameters or head observation wells (Yeh WW-G. Wat Resour Res 1986;22:95–108).All these facts demonstrate the potential of the algorithm to solve inverse problems of more complicated non-linear aquifer models naturally and quickly on the basis of finding suitable forms of the scaling matrix and the preconditioner.     

A viscous boundary for transient analyses of saturated porous media

An absorbing boundary for saturated porous media is developed that can be used for transient analyses in the time domain. The elastic constitutive equations for the saturated porous media follow Bowen's formulation. The method consists of applying viscous tractions along the artificial boundary. The absorbing boundary behaviour is assumed linear and isotropic. Hadamard's conditions provide the speeds of the dilatational and shear waves that propagate in saturated porous media. Since these expressions are frequency independent, the intensities of the viscous tractions are evaluated in the time domain, and the two dilatational waves are accounted for. The viscous tractions are defined from the drained characteristics, assuming an infinite permeability, at variance with the traditional ‘undrained’ method based on undrained characteristics and a null permeability. Solid media and materials with low permeability are also retrieved as subcases. The results show that, at no additional cost, this ‘drained’ method is more accurate for all permeabilities than the ‘undrained’ method, which disregards the existence of the second dilatational wave. Copyright © 2003 John Wiley & Sons, Ltd.     

Reduction in ground vibrations by using shaped landscapes

Reduction in traffic-induced ground vibrations by the use of shaped landscapes is investigated here by shaping the landscape surrounding a high-tech facility, using the landscape thus produced as a wave obstacle. The effects of the geometric parameters of a shaped landscape were examined in parametric studies. An architectural landscape design was also investigated in terms of its effectiveness in reducing traffic-induced ground vibrations. Finite element models, analysed in the frequency domain, were employed. The models involve a layer of soil and the underlying bedrock. It was found that anywhere from an appreciable reduction to an appreciable amplification of the vibrations produced can occur, depending upon the geometric parameters of the shaped landscape involved. The most effective shape was found for a topography that acted as a waveguide that reduced the level of vibration by approximately 35%.     

中国东部地区地震活动与构造应力场关系的有限元数值模拟

运用有限元数值模拟方法结合板块构造边界条件和板内构造分布特征对中国大陆东部地区的基本构造应力场进行数值模拟,根据研究区受周边板块作用的情况分析板块边界作用力的变化对内部应力场的影响。对比有限元数值模拟结果解释了中国大陆东部地区地震活动空间分布的特点。     

一种基于亥姆霍兹分解的大地电磁测深有限元正演预条件解法

本研究针对大地电磁测深法有限元数值模拟中,迭代法求解线性方程组效率较低的问题,利用亥姆霍兹分解原理,将电场矢量双旋度方程的预条件问题转化为基于矢量位的泊松问题和基于标量位的拉普拉斯问题,并在四面体非结构化棱边元离散的情况下,借助节点元辅助网格离散上述预条件问题,进一步利用代数多重网格方法（AMG）实施求解,最终实现预条件算法.利用经典的COMMEMI理论模型进行试算并与前人的积分方程解进行对比,验证了本文数值模拟程序与预条件方法的正确性和可靠性.此外,利用不同自由度规模的实验模型对这一预条件算法的效率进行了测试.结果表明,这一算法可以有效地提升大地电磁测深法棱边有限元数值模拟迭代法的收敛性,计算效率较通用的不完全LU分解预条件算法明显更高;在较大自由度网格（>1000万）数值模拟计算中,其算法效率及内存占用相对直接解法有较大优势,也使小型工作站上利用较大自由度的有限元网格进行大地电磁测深数值模拟计算成为可能.

     

径向分层TI孔隙介质井孔中激发的模式波的数值研究

横向各向同性介质是地层中普遍存在的一种各向异性介质.本文对径向分层TI孔隙介质包围井孔中激发的斯通利波和弯曲波的传播特性进行了理论计算,发现模式波在低频时更多的是反应原状地层的信息,而随着频率的增加侵入带参数逐渐起控制作用;Biot理论描述的地层衰减比速度更容易受井壁附近地层参数的影响.利用灵敏度曲线定量研究了不同频率下地层各个参数对相速度和衰减系数的贡献大小,主要结果显示模式波的衰减受水平渗透率影响明显,而垂直渗透率的变化对模式波几乎无影响;斯通利波对水平向传播的横波速度比弯曲波的灵敏度高.从单极子和偶极子声源在井孔中激发的全波波形也可发现,声波测井仪器较宽的声源频带和合适的源距设置有利于对不同径向深度上的地层声学参数进行成像.     

基于电场总场矢量有限元法的接地长导线源三维正演

电磁场数值模拟的背景场/异常场算法是三维正演的有效策略之一,优点为采用解析法计算电磁场背景场代替场源项、克服了场源奇异性,缺点为不适用于发射源布置于起伏地表或背景模型复杂的情形.总场算法是直接对电磁场总场开展数值模拟,其难点是有效加载场源、保证近区与过渡区数值解精度.本文以水平电偶源形式分段加载接地长导线源,并以电场总场Helmholtz方程为矢量有限元法控制方程,实现了基于非结构化四面体网格剖分的接地长导线源频率域电磁法三维正演.通过与均匀全空间中水平电偶源产生的电场解析解对比,验证了本文算法的正确性,并分析了四面体外接圆半径与其最短棱边的最大比值和四面体二面角最小值对数值解精度的影响规律.通过与块状高导体地电模型的积分方程法、有限体积法和基于磁矢量势Helmholtz方程的有限元法数值解对比,进一步验证了本文算法正确性,同时说明了非结构化四面体网格能够更加精细地剖分电性异常体,利于获得精确数值解.     

多个地震台网相互融合测定震级的贝叶斯算法

流动地震台网由于台站台基响应、震级测定方法和操作方式等多种因素的影响, 测定震级与区域内的固定地震台网震级测量结果存在差异, 从而使其记录到的而标准台网没有记录到的大量局域小地震不能直接补充到标准目录中.即使流动地震台网的震级测定公式中的量规函数经过校正, 两者在具体地震事件的震级测定上的差异仍不可避免, 因此急需解决这一普遍性问题的快速计算方法.本文先假定震级的观测值服从正态分布, 其中对固定台网而言, 测定震级的平均值为未知的真实震级, 对流动台网而言, 测定震级的平均值为真实震级的一个线性函数, 然后基于贝叶斯原理, 推导得出了不同台网观测震级的后验概率密度函数.应用于实际观测资料时, 可用最大似然估计来获得该算法中的矫正参数, 并计算后验的地震震级, 即矫正震级.我们将该算法实际应用于西昌流动地震台网(XC)与中国地震台网(CN)的震级融合问题中.基于矫正震级与原始测定震级的对比研究表明, 矫正震级与固定地震台网原始测定震级之间存在一定偏离度的线性对应分布, 而西昌流动地震台网测定震级相比矫正震级整体偏小.本文所提出的多台网震级融合的贝叶斯算法为流动台网测定的不精确震级提供了一种震级矫正的快速算法, 考虑固定台网与流动台网测定震级的一致性, 固定台网测定震级的稳定性和准确性也得到同步提高.

     

基于卷积神经网络和MPGA-LM算法的阵列侧向测井快速反演方法

在大斜度井/水平井环境下,阵列侧向测井受钻井液侵入、地层倾角和各向异性等多种因素影响,导致测井响应复杂,需借助反演手段提取地层真实电阻率.然而,阵列侧向测井三维正演效率低,难以满足测井资料快速反演和油气藏快速评价的需求.为此,本文基于深度学习并联合混合多种群遗传(MPGA)与列文伯格马奎特(LM)算法建立了一种快速反演...     

Describing near surface, transient flow processes in unconfined aquifers below irrigated lands: A deforming finite element model for heterogeneous aquifers

Hydrologic models of irrigated lands generally adopt either a basin-scale or a root-zone perspective. While basin-wide macro-scale models rely on the aggregation of important spatial and temporal data across large areas, micro-scale root-zone models depend on the definition of rigid boundaries around the zone of plant–soil–water interaction. In reality, irrigation management decisions are made on a field by field basis and can interact across field boundaries. This paper first describes a shallow water table model, based on deforming finite element (DFE) framework, to characterize the near-surface field-to-field hydrologic response to various irrigation and drainage management regimes along a gently sloping alluvial fan. The model is then enhanced through changing geometry of a fluctuating water table below a series of irrigated fields. Such an enhancement also offers computational flexibility relative to the saturated–unsaturated models commonly used in micro-scale studies. The model is designed with the alluvial fan aquifers of California’s western San Joaquin Valley as reference systems.     

A FEniCS-based programming framework for modeling turbulent flow by the Reynolds-averaged Navier-Stokes equations

Finding an appropriate turbulence model for a given flow case usually calls for extensive experimentation with both models and numerical solution methods. This work presents the design and implementation of a flexible, programmable software framework for assisting with numerical experiments in computational turbulence. The framework targets Reynolds-averaged Navier-Stokes models, discretized by finite element methods. The novel implementation makes use of Python and the FEniCS package, the combination of which leads to compact and reusable code, where model- and solver-specific code resemble closely the mathematical formulation of equations and algorithms. The presented ideas and programming techniques are also applicable to other fields that involve systems of nonlinear partial differential equations. We demonstrate the framework in two applications and investigate the impact of various linearizations on the convergence properties of nonlinear solvers for a Reynolds-averaged Navier-Stokes model.     

A quadratic element method for evaluating groundwater flow by the inversion of surface tilt with application to the Tono Area, Japan

By modifying a previous method with constant elements, we developed a quadratic element method for more accurately estimating groundwater flow by the inversion of tilt data. In this method: (1) a region of groundwater flow is divided into quadratic elements in which the change in groundwater volume per unit volume of rock (Δv) and the Skempton coefficient (B) vary in a quadratic manner with the coordinates, (2) the values of Δv are set to zero at the boundaries of the region of groundwater flow and (3) the sum of the squared second derivatives of Δv is adopted as a constraining condition that is weighted and added to the sum of the squared errors in tilt. First, analyses were performed for a flow model to determine the accuracy of this method for estimating groundwater flow and also to clarify the effect of the assumed size of a region of groundwater flow. These analyses showed that the quadratic element method proposed in this study gives a much better estimation of Δv than the constant element method and that a large region of groundwater flow should be assumed, rather than a small region, since the values of Δv at points outside of the actual region of groundwater flow are estimated to be nearly zero when a large region is assumed while these values are greatly overestimated when an excessively small region is assumed. Finally, the quadratic element method was applied to the site of the Mizunami Underground Research Laboratory in the Tono area, Japan. Inverse analyses were performed for tilt data measured by four tiltmeters with a resolution of 10⁻⁹ radians during the excavation of two shafts under the assumption that the rock mass is an isotropic and homogeneous half- space. The results showed that the method proposed in this study reproduced the tilt data very accurately. Thus, the distribution of Δv was estimated without sacrificing the reproducibility of the tilt data. The contour maps of B(1 + ν)Δv (ν: Poisson’s ratio) showed that the heterogeneous flow of groundwater occurred at the site and that groundwater volume decreased mainly in the area surrounded by two faults. The latter result is consistent with the finding obtained by previous investigations that these faults have low permeability in the direction perpendicular to the strike and may act as a flow barrier.