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

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

     

Three-dimensional finite-element modelling of magnetotelluric data with a divergence correction

A finite-element method for computing the electric field in a 3-D conductivity model of the Earth for plane wave sources, thus enabling magnetotelluric responses to be calculated, is presented. The method incorporates in the iterative solution of the electric-field system of equations the divergence correction technique introduced for finite-difference solutions by Smith (1996). The correction technique accelerates the development of the discontinuity of the normal component of the approximate electric field across conductivity discontinuities. The convergence rate of the iterative solution is improved significantly, especially for low frequencies. The correction technique involves computing the divergence of the current density for the approximate electric field, computing the static potential whose source is this divergence of the current density, and ‘correcting’ the approximate electric field by subtracting from it the gradient of the potential. This is repeated at regular intervals during the iterative solution of the electric-field system of equations. For the method presented here, the Earth model is discretised using a rectilinear mesh comprising uniform cells. Edge-element basis functions are used to approximate the electric field and nodal basis functions are used to approximate the correction potential. The Galerkin method is used to derive the systems of equations for the approximate electric field and correction potential from the respective differential equations. A bi-conjugate gradient solver was found to be adequate for the system of equations for the correction potential; a generalised minimum residual solver was found to be better for the electric-field system of equations. The method is illustrated using the COMMEMI 3D-1A and 3D-2A models.     

三角网格有限元法波动模拟的数值频散及稳定性研究

三角网格有限元法能够准确模拟复杂构造和复杂介质条件下的地震波场,数值频散和稳定性条件是地震波数值模拟中参数选择的主要依据.基于均匀的线性三角网格单元,根据结构刚度矩阵的组装原理以及平面波理论,推导了集中质量矩阵下两种网格结构的声波频散函数以及稳定性条件,并对数值频散特性以及稳定性进行了详细研究:三角网格单元中波动的数值频散除了受到空间采样间隔、单元网格纵横比和波传播方向等常规因素的影响外,还受到网格布局的影响,过锐或过钝的三角单元会对波动数值频散产生不良的影响,不同类型的单元网格、单元纵横比对应着不同的稳定性条件,正三角单元中的波动具有较好的数值频散特性,其数值各向异性(频散随波传播方向的变化)效应最弱,稳定性条件也较为宽松.最后通过数值模拟直观地验证了以上分析结果,为有限元正演三角网格的剖分和参数的设置提供一定的理论依据.     

基于自适应有限元正演的大地电磁法三维反演算法研究

目前大地电磁法三维反演大多基于有限差分正演.反演使用规则六面体网格无法有效模拟复杂地形, 同时正反演网格同套存在严重影响反演可靠性的问题.针对上述两个问题, 本文利用自适应有限元算法, 开发了基于非规则六面体的高精度大地电磁三维正演方法; 反演中, 使用独立的反演网格和正演网格来提高正演计算精度和反演可靠性.这一思路既保证了正演响应和灵敏度矩阵计算的精确性, 又降低了因反演参数过多造成的不唯一性.最后, 通过算例验证了正演算法的精确性和反演算法对地形处理的有效性.

     

基于矢量有限元的带地形大地电磁三维反演研究

研究了基于矢量有限元方法的大地电磁带地形三维反演算法并开发了三维反演计算程序代码.在大地电磁场正演数值模拟方面,采用并行直接稀疏求解器PARDISO且无需进行散度校正的快速正演方案,对典型地形模型,在中等规模计算条件下,与双共轭梯度法（BICG）计算结果比较,发现PARDISO比BICG快10倍以上;通过理论模型试算,并与前人的有限元法计算结果对比,验证了带地形三维正演计算程序的正确性.在反演方面,本研究基于共轭梯度方法编写了大地电磁带地形三维反演代码,为了避免直接求取雅可比矩阵,将反演中的雅可比矩阵计算问题转为求解两次"拟正演"问题,进而将PARDISO的快速正演方案应用于"拟正演"问题的求解,以提高反演计算效率.利用开发的反演算法对多个带地形地电模型的合成数据进行了三维反演,反演结果能很好地重现理论模型的电性结构,验证了本文开发的三维反演算法的正确性和可靠性.最后,利用该算法反演了某矿区大地电磁实测数据,反演得到的三维电性结构清晰地反映了研究区的地电特征,将反演结果与该区已有地质资料结合进行解释,应用效果明显,进一步验证了本文算法的有效性.     

基于二次场方法的并行三维大地电磁正反演研究

快速且高精度的三维大地电磁法正反演是目前研究的热点.由于大地电磁法场源的平面波特性,以往的正演方法大多采用直接求解总场的方法,在边界强加二维边界条件.本文提出了一种基于二次场方法的三维大地电磁法正演算法,将平面波在层状背景模型中的响应作为场源项,得到二次场满足的偏微分方程,并利用交错网格有限差分法求取二次场.与其他学者的基于总场方法的结果的对比证明了本文采用方法的正确性.在基于二次场的正演算法基础上,实现了基于L-BFGS的三维反演方法,并对公开的数据集进行了反演.另外,针对大地电磁法的多频率观测特性,采用了基于MPI的分频并行策略对程序进行并行化,可达到接近线性的加速比.

     

大地电磁三维交错网格有限差分数值模拟中的散度校正方法研究

散度校正对加快大地电磁三维正演速度、提高计算精度均具有重要意义.本文基于交错网格有限差分开展了三维大地电磁正演模拟,从Maxwell方程中磁场散度为零的特征出发,提出了三种散度校正方法:直接磁场散度校正(DDCM)、磁场散度残差复校正(RPCM)、磁场散度残差实校正(RRCM),并将上述校正方法用于垂直接触带、典型低阻异常体和国际标准地电模型的正演计算.结果表明:相对于其他两种散度校正和不做散度校正求解方案,磁场散度残差复校正具有计算时间短,收敛速度更快,计算结果更精确的优点.此外,不同网格剖分下的磁场散度残差实校正方案计算结果显示出该校正方案具有良好的稳定性和适用性.     

最优系数有限单元法色散介质GPR频率域正演模拟

工程实际勘探对象如土壤、岩石等多为色散介质,雷达波在其中传播时易发生衰减与畸变,应用常规有限单元法（Finite Element Method,FEM）方法进行数值模拟时,存在数值频散现象.为此,作者以色散介质为研究对象,开展最优系数有限单元法探地雷达（Ground Penetrating Radar,GPR）频率域正演.首先,分析了有限元质量、刚度矩阵的约束条件对有限元求解精度的影响,基于归一化相速度与1的误差最小策略,利用最小二乘法,仅需三个优化参数求取最优的有限元刚度矩阵与质量矩阵.四种不同方法的频散曲线分析及精度对比实验结果表明,优化矩阵在单位波长仅需4.8个网格点下便可达到误差小于0.2%的精度;而一致、集中和折衷矩阵不仅需要更多的网格点,且误差较大.然后,将精确完全匹配层（Exact Perfectly Matched Layer,EPML）吸收边界条件引入最优系数频域有限单元（Finite Element Frequency Domain,FEFD）算法中,简化了吸收参数优化过程,取5层即可达到常规完全匹配层（Perfectly Matched Layer,PML）的10层的吸收效果,能够有效提升正演效率.并将基于EPML的最优系数有限单元法算法引入到城市道路病害模型正演中,实验表明：本文算法能有效压制频散并实现实际色散介质高精度模拟,模拟结果更接近波在地下介质中的实际传播特性.

     

基于第三类边界条件的大地电磁三维矢量有限元任意各向异性正演

地球内部介质的各向异性对地球物理场解译有很大影响, 研究各向异性介质中大地电磁响应具有重要的意义.边界条件是影响电磁场正演精度的一个关键因素, 其中第一类边界条件需要将底面边界设置在离异常体足够远的地方, 面临着计算规模大、求解速度慢的问题.相比第一类边界条件, 第三类边界条件具有计算规模更小、计算精度更高的优点, 在三维各向同性正演中被广泛使用.然而, 各向异性使得第三类边界条件理论变得更为复杂, 目前尚未见到采用第三类边界条件的大地电磁三维各向异性正演.本文推导出各向异性介质中三维矢量有限元底界面的第三类边界条件, 通过模型计算对算法的可靠性和精度进行了验证.结果表明, 当精度相同时, 第三类边界条件的计算规模可以更小; 而当底界面离异常体较近时, 第三类边界条件有着更高的精度.

     

耦合PML边界条件的三维大地电磁二次场有限差分法

本文将大地电磁场分解为一次场和二次场,应用交错网格有限差分法模拟计算大地电磁二次场,并引入各向异性最佳匹配层（PML）吸收边界条件作为二次场边界条件,实现了耦合PML吸收边界条件的三维大地电磁二次场有限差分正演模拟.为了确保正演的稳定性和效率,QMR求解器和磁感应矢量散度校正技术被用于PML吸收边界条件下系数矩阵的快速求解.三维模型正演响应表明,基于二次场的三维大地电磁有限差分算法具有较高的计算精度和可靠性.通过计算分析不同PML吸收因子条件的大地电磁正演结果,显示在适当的吸收因子下,PML吸收边界条件可较大幅度的减小外边界距离,从而有效的压缩模型求解空间,最终提高三维大地电磁正演模拟的效率.

     

有限元感应测井模拟的背景场选择方法研究

在感应测井的有限元模拟中,为了消除源的奇异性,一般将总场分解为背景场和散射场.本文定量研究了不同背景场选择方法对计算精度的影响.首先,分析了在均匀和径向分层介质中不同背景电导率对长短源距线圈系有限元模拟结果的影响;其次,利用三层介质模型对比了选择源点附近和线圈系中点附近地层电导率作为背景电导率的结果.研究结果表明,如果在选择背景场时,只考虑源附近的散射场梯度而不同时考虑源和接收器附近的散射场梯度,计算误差明显增大.在此基础上提出一种利用Gianzero几何因子计算的视电导率作为背景电导率的新的背景场方法.该方法综合考虑了围岩、井眼、线圈距等因素,特别是在介质分界面处,可有效减小计算误差,并取得了满意的精度.本研究为复杂环境下的感应测井模拟的背景场选择提供了指导和依据.     

修正辛格式有限元法的地震波场模拟

三角网格有限元法具有网格剖分的灵活性,能有效模拟地震波在复杂介质中的传播.但传统有限元法用于地震波场模拟时计算效率较低,消耗较大计算资源.本文采用改进的核矩阵存储（IKMS）策略以提高有限元法的计算效率,该方法不用组合总体刚度矩阵,且相比于常规有限元法节省成倍的内存.对于时间离散,将有限元离散后的地震波运动方程变换至Hamilton体系,在显式二阶辛Runge-Kutta-Nystr m（RKN）格式的基础之上加入额外空间离散算子构造修正辛差分格式,通过Taylor展开式得到具有四阶时间精度时间格式,且辛系数全为正数.本文从理论上分析了时空改进方法相比传统辛-有限元方法在频散压制、稳定性提升等方面的优势.数值算例进一步证实本方法具有内存消耗少、稳定性强和数值频散弱等优点.

     

基本结构有限元算法及大地电磁测深一维连续介质正演

提出了一种新的数值方法——基本结构有限元法. 从经典的伽辽金法（Galerkin method）出发,首先定义了基本结构插值基函数,在此基础上建立了基本结构方程,然后用有限元法进行进一步的具体实现. 该算法被成功应用到大地电磁测深一维连续介质正演计算中. 采用了6种不同的单元插值基函数进行计算,并对计算结果进行了比较和分析.     

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

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

龙滩水库诱发地震三维孔隙弹性有限元数值模拟

本文以龙滩水库为例,根据库区地质构造、深部速度结构及数字地面高程,建立了库区三维有限元模型,基于孔隙弹性理论计算了水库蓄水过程中库底断层和围岩体孔隙压力、有效附加正应力、剪应力和库伦应力的动态变化,并结合水库蓄水后库区地震活动时空分布的特征,讨论了RIS时空演化与库水加卸载及渗透过程的动态响应关系及其可能的成因机制.结果表明：（1） 龙滩水库蓄水后地震活动呈现出明显的丛集性,主要分布在罗妥（丛Ⅰ）、八茂（丛Ⅱ）、拉浪（丛Ⅲ）、坝首（丛Ⅳ）和布柳河（丛Ⅴ）5个水库蓄水后淹没的深水区,这些区域也恰恰是库水加卸载及渗透过程中ΔCFS增加最明显的区域,而ΔCFS的影区几乎没有地震发生,表明水库蓄水后库区地震活动与ΔCFS的变化密切相关.（2） 在水库蓄水过程中,与水库有直接水力联系且渗透性较好的断裂成为地表水体附加水头压力向深部扩散的优势通道,沿此通道附加水头压力扩散的最大深度达13 km左右,震旦系—古生界以碳酸盐岩为主的地层成为附加水头压力扩散的主体层位,这与蓄水后库区中、小地震震源深度均小于13 km,且优势分布在5~10 km的特征相吻合,表明由于孔隙压力的存在降低了岩石的抗剪强度,同时部分抵消了围压的影响,致使该层位的岩体易于产生脆性破坏从而诱发地震活动.（3） 无论是深部还是浅部,各丛地震密集发生的时段绝大部分与相应深度ΔCFS加速升高或阶段性高值时段相重叠,可能说明在库水位快速抬升或阶段性高值时段,受外部荷载加载速率快速升高的影响,库底岩体和断层、裂隙等结构面更容易实现失稳扩展;深、浅部地震响应时间、活动频度和强度的差异可能与不同层位岩体力学性质及渗透性能的不均匀性有关.（4） 各丛地震诱发的物理力学机制有所不同.丛Ⅰ、丛Ⅱ、丛Ⅲ地震的诱发可能与库体重力荷载、孔隙压力扩散和库水浸润弱化3种作用都有关;丛Ⅳ地震的诱发主要受控于库体重力荷载作用,孔隙压力扩散和库水浸润弱化不起主导作用;丛Ⅴ地震的诱发主要受孔隙压力扩散和库水浸润弱化作用的影响,库体重力荷载作用一定程度上抑制了地震的发生.     

自适应非结构有限元MT二维起伏地形正反演研究

在山区进行MT勘探时,用规则网格有限元方法模拟起伏地形会受到限制.本文采用非结构三角网格可以有效地模拟任意二维地质结构,如起伏地形、倾斜岩层和多尺度构造等.正演引入自适应有限元方法,其在网格剖分过程中能根据单元误差自动细化网格,保证了正演结果的精度.将自适应有限元与Occam算法结合,且引用并行处理技术提高正反演计算速度.通过对比两个理论模型,讨论了地形对MT正演响应的影响;其次进行了不同地电模型带地形反演展示了本文算法的正确性和适用性;最后将该方法应用于实测MT数据处理,证明了自适应非结构有限元方法是复杂地形下处理MT数据的有力工具.     

可控源电磁三维频率域有限元模拟

本文采用电磁场的磁矢量位和电标量势,将Maxwell方程组化为位势的类似于Helmholtz型方程,并引入罚项及稳定化方法克服了电磁三维有限元算法中的伪解及数值不稳定性;采用人工边界把计算区域局域化,将均匀半空间中水平电偶极子源产生的位势值作为人工边界上的第一类边界条件以表示源的作用,减少了实际的计算区域.理论模型和复杂模型的计算结果均表明,可控源电磁三维有限元数值模拟给出了稳定、可靠的电磁场分布.     

三角网格有限元法声波与弹性波模拟频散分析

本文对声波与弹性波方程进行有限元法离散,构造有限元法频散关系的一般特征值问题,分析了时间离散格式为中心差分的三角网格有限元法声波与弹性波模拟的频散特性. 比较了三种质量矩阵即分布式质量矩阵、集中质量矩阵和混合质量矩阵对有限元法频散的影响;选取四种典型三角网格,分析了混合质量矩阵有限元（MFEM）频散的方向各向异性;数值频散、方向各向异性随插值阶数的增加逐渐减弱,当空间为三阶插值时,频散主要表现为随采样率的变化而几乎无明显方向各向异性, 其频散幅值也较小. 控制其他影响因素不变的情况下,研究了不同波速比介质中弹性波的数值频散. 最后给出了三角网格MFEM的数值耗散性.     

复电阻率法二维有限元数值模拟

伴随着复电阻率法的广泛应用,发展精确和快速的正演和反演算法成为复电阻率法研究的重点.本文采用基于三角单元剖分的有限单元法进行了复电阻率二维数值模拟研究.为了提高计算速度,对无穷远边界进行了近似处理.整个正演计算过程分为两步,首先采用有限单元法计算四个不同频率的视复电阻率数据,然后对前一步得到的视复电阻率数据采用递推算法计算视Cole-Cole参数.采用这种正演算法与一维正演的结果进行了对比,验证了本文方法的正确性.设计了两个二维极化模型,数值模拟结果表明视复电阻率和Cole-Cole视参数等值线断面图对于异常目标体都有比较明显的反映.