基于非结构化网格的大规模磁测数据正则化反演研究.pdf

第 39 卷第 5 期2023年9月UraniumGeology铀矿地质Vol.39 No.5Sep.2023DOI:10.3969/j.issn.1000-0658.2023.39.073基金项目 国家自然科学基金项目（编号：42164008）和江西省自然科学基金项目（编号：2019BAB202012）联合资助。收稿日期 2023-03-02 改回日期 2023-05-19作者简介 刘阳（1996），男，甘肃陇西人，硕士研究生，研究方向磁法勘探。E-mail：通信作者 张志勇（1978），男，河北承德人，教授，博士，研究方向地球物理数值模拟与反演。E-mail：基于非结构化网格的大规模磁测数据正则化反演研究刘阳，张志勇（东华理工大学 地球物理与测控技术学院，江西 南昌 330013）摘要 磁测数据解析正演的内存需求与模型规模、观测点数成正比，以其为依托开展大规模磁测数据反演时，计算成本较大。文章利用非结构化网格有限单元法进行磁测正演，分析了地表观测网格、边界节点数和局部网格密度等因素对正演计算精度的影响；在此基础上，基于灵敏度矩阵隐式存储的高斯-牛顿法求解正则化反演目标函数，实现了大规模磁测数据三维反演。理论模型试算和铀矿区域实测数据案例表明，文章所述方法可保证正演计算的精确度，提高了计算效率、减少了反演的内存消耗，可用于复杂条件下大规模磁测数据的三维反演。关键词 磁法勘探；磁化率；有限单元法；正则化反演文章编号 1000-0658(2023)05-0831-12 中图分类号 P631.2 文献标志码 A磁法勘探已广泛应用于铀矿勘探、深部地质构造圈定、地热调查等领域1-4，大规模铀矿地质勘探亟需准确、高效的磁法解译。位场正演是磁法勘探解释任务的基础，其计算成本依旧是阻碍三维反演发展的一个严重问题5。位场正演模拟主要有空间域和频率域两类方法，其中空间域方法在处理起伏地表与非规则离散模型等方面具有较大优势。空间域方法又可分为解析法和数值法。通过闭合积分公式建模的空间域方法已有众多研究成果6-7，简单矩形棱柱积分公式难以拟合复杂地质体，均匀多面体剖分算法逐渐成为位场正演建模的主流8。然而，三维复杂地质模型计算公式复杂，适应性较差、计算量较大9，以解析法正演开展大规模离散反演时需要大量积分运算和内存消耗，计算成本巨大，因此，位场解析法建模的应用受到了限制。位场数值法建模主要有有限差分法10、有限体积法11-12、积分方程法13、有限单元法等14，与解析法相比，数值法计算精度略低，但计算时间和内存消耗方面具有一定优势。其中有限差分法通常采用结构化网格，但在物性参数复杂分布或场域分布几何特征不规则时适应性差，存在一定局限性。有限单元法既可采用结构化网格，也可采用非结构化网格。与结构化网格相比，非结构化网格在要求较高精度或可接受较低精度的区域，网格的局部细化或粗化可以在不影响这些区域以外的网格的情况下执行12。因此，使用非结构化网格的数值格式会有更高的精度、更低的计算时间和内存占用。有限单元法以其适用性强、可以模拟任意复杂地形及地质体而受到广泛运用15-17，刘阳铀 矿 地 质第 39 卷如地震18、直流电阻率15、自然电位19、大地电磁20正反演中均得到应用。前人基于有限单元法开展了大量位场正反演的研究，Zhang等21将遗传算法与有限单元法相结合实现了一种替代方法反演重力异常数据，重建地下三维密度结构；Cai 等14通过一种快速的有限单元法来解决复杂地质条件下引力位的边值问题；Maag 等22由重力位的偏微分方程出发推导 出 基 于 有 限 单 元 法 的 正 演 和 反 演 问 题；Jahandari 和 Farquharson12采用有限单元法和有限体积法与解析法对比实现了重力位正反演；但当前文献少有以非结构化有限单元法开展磁场正反演的研究。大规模磁测数据反演中，灵敏度矩阵的储存与计算一定程度上制约着计算规模的大小。通过将灵敏度矩阵及其转置与向量乘积转换为正演过程，进行有限内存最优化求解的反演策略已应用到地球物理反演问题中，其中算法包括共轭梯度法23、有限内存拟牛顿法24-25、高斯-牛顿法26等，与其他最优化算法相比，高斯-牛顿法在保证快速收敛的同时，有效减少了计算量，因此，高斯-牛顿法适用于大规模磁测数据反演。针对以解析法正演为依托的磁异常反演及规则网格数值解对复杂地质体拟合程度不足、计算成本大等缺点，本文采用非结构化网格有限单元法进行磁异常正演，通过灵敏度矩阵隐式存储的高斯-牛顿方法求解反演目标函数实现磁测数据正则化反演。研究工作首先采用局部网格加密技术和最小二乘空间梯度计算方法提高有限单元正演计算精度，通过与解析法对比，对非结构化网格有限单元法正演精度的影响因素进行了试算讨论；在此基础之上，开展了理论模型与铀矿区域实测数据的反演，验证了本文提出的三维磁测数据反演算法的可行性。1三维有限单元正演1.1正演理论在假设无传导电流的空间域，磁位基本方程可表示为：(0(-U-U)+Mr)=0（1）式中：0为真空磁导，H/m；U为磁位；为磁化率；为哈密尔顿算子；Mr为剩余磁化强度矢量，A/m。当在两种介质存在的空间域中，采用有限单元法求解方程式（1），考虑到两种介质分界面磁位连续且无散场的特点，得到公式：U1=U2（2）1U1n+2U2n=Mr1，n+Mr2，n（3）式中：n为分界面上的单位正法线方向；1和2为边界两种介质中的磁导率；U1和U2为边界两种介质中磁位；Mr1和Mr2为边界两侧介质中剩余磁化强度矢量，Mr1，n和Mr2，n是Mr1和Mr2在外法向的投影，则有第一类边界条件：U=-T r（4）式中：T为总磁场强度矢量，T；r为坐标原点到无穷远边界距离矢量。第二类边界条件可表示为：Un=-Tn0（5）磁场边值问题与求解下式变分问题等价：F(U)=12(U)2-0U Mr d+TnUdF(U)=0（6）式中：Tn是磁场T在边界外法线方向上的投影，为研究区域。对计算区域进行非结构化四面体离散，在每个剖分单元采用线性函数插值，对各单元泛函求和，得到下列线性系统：KU=P（7）式中：K为刚度矩阵（对称稀疏矩阵）；U为节点磁位；P为边界条件形成源项。求解联立方程组可得到各节点磁位U，通过数值微商可得到各磁场分量，但由于非结构化网格中网格大小不一致，基于线性插值计算磁位，对磁位进行求导磁场为常数，使得在计算磁场时的精度和稳定性存在问题。本文拟采用与观测点共点的所有四面体进行空间梯度的最小二乘计算。假设定义观测点坐标为(x0，y0，z0)，相应磁位为U，与观测点 共 点 的 所 有 四 面 体 n 个 顶 点 坐 标 定 义 为(xi，yi，zi)，相应磁位为Ui，其中i=1，2n，则有：832刘阳，等：基于非结构化网格的大规模磁测数据正则化反演研究第 5期U(x，y，z)=ax+by+cz+f（8）节点组成方程组为：x0y0z01x1y1z11xnynzn1 abcf=U0U1Un（9）可将式（9）表示为：AL=U（10）令E=(U-AL)T(U-AL)，对L进行求导可得：EL=(U-AL)T(U-AL)L=2AT(AL-U)（11）令上式为 0，得L=(ATA)-1ATU，求解方程可得到坐标点处梯度为：U=(x0，y0，z0)=ax+by+cz（12）进而可得到观测点处磁场强度。1.2正演精度影响因素分析正演精度直接影响反演的可靠性，对于非结构化有限单元正演，其精度受多种因素影响。本文从观测网格密度、边界节点数两方面进行讨论，并与均匀多面体的闭合积分公式27计算出的总磁异常进行精度对比。设计模型如图 1所示，假设在地下空间磁化率为 0.006 SI的背景中，设置一个磁化率为 0.1 SI的异常体，异常体顶板埋深为 500 m，结构为 1 000 m1 000 m1 000 m，地磁场强度在 X、Y、Z三个方向的分量分别为（30 000 nT，0 nT，40 000 nT）。选取观测区任意部分区域测点验证有限单元法的计算精度，计算区域大小为 250 m250 m。图 1 磁化率模型Fig.1 Magnetic susceptibility modela俯视图；b立体图。1.2.1观测网格影响在有限单元计算中，以观测点作为剖分网格中心点进行剖分，因此，有限单元计算精度依赖于剖分网格到中心点的距离。对不同观测网格产生总磁异常进行计算，边界节点数为 38，计算结果如图 2 所示。不同观测网格的剖分网格属性及各类误差信息如表 1 所示。随着观测网格变密，平均绝对误差从 3.67 nT降至 1.60 nT，最大绝对误差从 14.08 nT 下降至 5.56 nT，平 均 相 对 误 差 从 0.047 下 降 至0.018，计算总磁异常精度明显提高。1.2.2边界节点数影响为分析边界节点数对有限单元正演模拟计算的影响，选取边界节点 15，25，35，45进行试算，观测网格为 50 m50 m，以 1.2 倍指数扩边得到图 3和表 2。计算结果显示，随着边界节点数的增加，平均绝对误差从 80.15 nT 下降至 3.77 nT，最大绝对误差从 95.33 nT 下降至 12.31 nT，平 均 相 对 误 差 从 0.460 下 降 至0.046，最大相对误差减小，计算精度明显提高，而过少边界节点计算总磁异常存在较大的误差。1.2.3局部网格密度影响观测网格密度和边界节点数对于有限单元计算精度均有影响，而观测网格变密导致计算效率的下降，此外，观测网格密度达到一 833铀 矿 地 质第 39 卷定程度，继续加密对于计算精度有较小的影响15。因此，本文采用网格加密的方式提高磁异常计算精度，网格加密的方法主要有局部 体 积 加 密 方 法（Local Volume-Refinement，简 称 LVR）、全 局 体 积 加 密 方 法（GlobalVolume-Refinement，简称 GVR）和局部测点加密 方 法（Local Node-Refinement，简 称 LNR）。其中 GVR 法加密基于结构化网格会导致过多的冗余加密节点，精度提高不明显；LVR 法是在加密需求区域无应有的加密，反而在其他区域加入过多节点，导致计算量增大；LNR法相对于 LVR 和 GVR 法可以减少内存消耗和计算时间，因此，本文采用 LNR 网格加密方式。局部测点网格细化策略是基于观测点位置下方插入节点，观测网格为 50 m50 m，扩图 2 不同观测网格下的总磁异常等值线及绝对误差图Fig.2 Contour of total magnetic anomaly and absolute error under different observation gridsa观测网格 50 m50 m等值线；b观测网格 25 m25 m等值线；c观测网格 12.5 m12.5 m等值线；d观测网格 50 m50 m绝对误差；e观测网格 25 m25 m绝对误差；f观测网格 12.5 m12.5 m绝对误差。表 1 地表不同观测密度所产生的网格参数及各类误差信息表Table 1 Grid parameters and error information of different observing density属性节点数四面体单元数三角单元数边数最大绝对误差平均绝对误差最大相对误差平均相对误差50 m50 m400 9442 552 7135 106 6292 954 85914.08 nT3.67 nT0.3720.04725 m25 m414 6722 637 6825 276 7353 053 72411.94 nT2.66 nT0.2470.03012.5 m12.5 m444 6352 823 0575 647 8493 269 4265.56 nT1.60 nT0.1540.018 834刘阳，等：基于非结构化网格的大规模磁测数据正则化反演研究第 5期边节点数为 45，得到计算结果如图 4 和表 3 所示。局部网格细化可提高正演计算精度，插入节点深度较浅时能获得更高精度的异常分辨率，当插入节点后剖分网格数量增加明显，而插入节点接近观测点时剖分网格数量明显多于插入节点较深的情况，当在观测点下 60 m、20 m、1 m 插入节点时平均绝对误差从 2.47 nT减小至 1.11 nT，平均相对误差从 0.028 下降至0.011，相比同一观测密度计算，局部网格加密后磁异常精度和稳定性明显提高。以上试算结果表明，通过提高观测网格密度、增多边界节点数及采用局部网格细化技术等可实现非结构化网格有限单元法计算磁异常精确解。1.3有限单元计算效率分析为说明非结构化网格有限单元法计算的低耗存和高效率的特点，本文测试平台为 Intel（R）Core（TM）i7-9700F CPU3.00 GHz，RAM32.0 GB。在同样的观测网格及相同剖分网格参数下，进行有限单元法与解析法计算内存和时长的统计（表 4 和图 5）。如表 4 和图 5 所示，与解析算法相比，同样的计算条件和剖分网格参数下，有限单元法计算时长与内存需求明显优于解析算法。图 3 不同边界节点数计算总磁异常等值线及绝对误差图Fig.3 Contour of total magnetic anomaly and absolute error calculated with different boundary nodesa边界节点数 15等值线；b边界节点数 25等值线；c边界节点数 35等值线；d边界节点数 45等值线；e边界节点数 15绝对误差；f边界节点数 25绝对误差；g边界节点数 35绝对误差；h边界节点数 45绝对误差。表 2 不同边界节点数计算产生网格参数及各类误差信息表Table 2 Grid parameters and error information calculated with different boundary nodes属性节点数四面体单元数三角单元数边数最大绝对误差平均绝对误差最大相对误差平均相对误差15363 7492 322 4744 645 6312 686 90595.33 nT80.15 nT0.7020.46025376 3582 399 3464 799 7102 776 72129.29 nT19.77 nT0.4300.17135397 9332 535 8195 073 0092 935 12212.41 nT4.32 nT0.2100.05045428 8052 727 3585 456 4703 157 91612.31 nT3.77 nT0.1960.046 835铀 矿 地 质第 39 卷2大规模磁测数据反演算法2.1正则化目标函数正 则 化 反 演 目 标 函 数 一 般 形 式 可 表示为28：(m，dobs)=d(m，dobs)+m(m)=Wddobs-Am22+Wm()m-mapr22（13）式中：d(m)为数据误差函数；m(m)为模型稳定函数，本次研究工作采用最小结构稳定函数29；为正则化因子，是平衡数据误差函数和模型稳定函数之间的权重系数，采用 L-Curve 方法选取30；dobs为观测数据向量；m为离散模型向量；mapr为先验模型参数向量；A为正演算子；Wd为数据协方差矩阵；Wm为模型权重矩阵。2.2目标函数优化本次研究工作采用高斯-牛顿方法求解目标函数：m(n+1)=m(n)+m(n)（14）式（13）数据误差项和模型稳定项展开可得：n+1d(m)=Wd(dobs-Am(n)-Am(n)22（15）n+1m(m)=Wm(m(n)+m(n)-mapr)22（16）将上述方程对m(n)求导，可得高斯-牛顿方程：A(n)TWTdWd(dobs-Am(n)-WTmWm(m(n)-mapr)=(A(n)TWTdWdA(n)+WmTWm)m(n)（17）图 4 不同深度插入节点计算磁异常等值线及绝对误差图Fig.4 Contour of magnetic anomaly and absolute error calculated by inserting nodes at different depthsa插入节点深度 60 m等值线；b插入节点深度 20 m等值线；c插入节点深度 1 m等值线；d解析公式计算总磁异常等值线；e插入节点深度 60 m绝对误差；f插入节点深度 20 m绝对误差；g插入节点深度 1 m绝对误差。表 3 不同深度插入节点产生网格参数及各类误差信息表Table 3 Mesh parameters and various error information calculated by inserting nodes at different depths属性节点数四面体单元数三角单元数边数最大绝对误差平均绝对误差最大相对误差平均相对误差60 m465 4142 952 1455 906 0273 419 2958.08 nT2.47 nT0.2010.02820 m467 0172 963 2075 928 1693 431 9784.74 nT1.50 nT0.0680.0141 m472 6542 998 9365 999 5993 473 3163.092 nT1.11 nT0.0280.011 836刘阳，等：基于非结构化网格的大规模磁测数据正则化反演研究第 5期式 中：WmTWm=sDTD+xWxTDTDWx+yWyTDTDWy+zWzTDTDWz，其中D为深度加权函数产生的对角矩阵，Wx、Wy、Wz分别为X、Y 和 Z 三个方向的粗糙度矩阵31。采用迭代法29求解式（17）。三维磁测反演中直接储存灵敏度矩阵占用大量内存消耗，因此，本文反演过程中采用灵敏度矩阵隐式储存的策略减少内存消耗。如下为实现过程，采用互换定理求取灵敏度，观测数据可表示为：H=U（18）采用J表示灵敏度矩阵，J中任意元素表示为：Jij=Himj=()Umji（19）式中：Hi表示第 i个观测数据；m表示反演模型参数；mj为第 j 个反演模型参数。将其与有限单元线性方程式（7）结合，则：KmU+UmK=Pm（20）Um=K-1()Pm-KmU（21）将上式代入式（19），可得：J=K-1()Pm-KmU（22）令 磁 位 空 间 梯 度 变 换 矩 阵 为G，S=()Pm-KmU，则灵敏度矩阵及其转置可表示为：J=GUm=GK-1S（23）JT=GTK-1ST（24）将灵敏度矩阵及其转置与任意矩阵的乘积转换成线性方程的求解，并且该方程与有限单元具有相同的矩阵：JV=GK-1SV（25）JTV=GTK-1STV（26）由式（25）和（26）所示，求取灵敏度矩阵（或转置）和任意矩阵的乘积可包含在灵敏度求解过程中，因此，实现了灵敏度矩阵的隐式表 4 有限单元法与解析法计算效率对比Table 4 Comparison of computational efficiency between finite element method and analytical method边界节点数202530354045解析解内存/MB3 0725 1227 16614 23818 63633 792时长/s1 008.8911 491.862 741.3754 464.3287 016.39110 330.157有限单元法内存/MB788395108125146时长/s7.2513.32323.20341.28368.625122.318节点数36 0214 4786 0327 6429 78312 167四面体单元数14 29314 70119 62323 76029 64836 561三角单元数31 08833 00844 14853 92367 39883 124边单元数20 39622 76930 56837 81147 43458 622图 5 计算内存与时长曲线对比图Fig.5 Contrast of calculation memory versus duration curve between the method of finite element and analysisa解析法计算时长与内存曲线；b有限元法计算时长与内存曲线。837铀 矿 地 质第 39 卷储存，减少了内存消耗，提高了计算效率。3模型试算分析3.1单一模型首先采用单一模型验证本文提出算法的可行性，单一模型和复杂模型正演数据中均加入 5%随机噪声。模型设置如图 1，采用总磁异常数据进行反演，观测网格为 50 m50 m，测点数为 4 900 个，初始正则化因子为 1.010-5，迭代 75 次。如图 6 所示，采用本文算法对单一模型反演取得了良好的效果，其中白色矩形框为模型位置，反演结果磁化率分布与真实模型位置有较好对应。从 RMS误差曲线可知（图7），随着迭代的进行，RMS 总体趋势稳定下降趋近于 1.0，说明反演稳定收敛；反演中通过将灵敏度矩阵及其转置与任意向量的乘积转换成正演计算，内存需求减少，计算效率明显提高。3.2复杂模型为进一步测试本文提出的反演方法对复杂模型的反演效果，设置复杂磁化率模型如图8 所示。假设在地下空间磁化率为 0.006 SI 的背景中，设置 4 个磁化率均为 0.1 SI 的矩形异常体，异常体结构及属性见表 5，观测网格为50 m50 m，测点数共计 4 900 个，地磁场强度与单一模型等同。对于复杂模型进行磁化率反演取得同样效果（图 9），反演磁化率分布位置与设定的 4个异常体模型位置均有对应，且反图 6 单一模型反演结果Fig.6 Inversion results of single modela深度 Z=1 000 m横切片图；bX和 Y方向纵切片图。演磁化率横向分辨率有较好体现（图 9a）。从复杂模型反演 RMS 误差曲线可以发现（图 10），RMS 随迭代次数逐渐下降并最终趋近于 1.0，说明反演稳定收敛；相比解析法，以有限单元法开展复杂模型离散反演计算效率明显提高，试算结果表明，本文提出反演算法对于复杂模型同样具有可行性。表 5 异常体属性Table 5 Attributes of anomaly body异常体1234磁化率/SI0.10.10.10.1上顶埋深/m100150200250边长/m200200300300400400500500图 7 单一模型反演 RMS误差曲线Fig.7 RMS error curve of inversion result by single model 838刘阳，等：基于非结构化网格的大规模磁测数据正则化反演研究第 5期图 10 复杂模型反演 RMS误差曲线Fig.10 RMS error curve of inversion result by complex model3.3实测数据反演磁法在铀矿地质勘查中应用广泛，为验证本文提出算法处理铀矿区域采集实际数据的有效性，选取某铀矿研究区带地形磁测数据进行反演。研究区岩性以花岗岩为主，勘探目标为辉绿岩，观测网格为 25 m25 m，测点数共计 3 025 个，研究区反演最高处为 15 m，磁场强度为当地地磁场，磁异常为辉绿岩引起。初始正则化因子为 1.010-8，共计迭代 28 次。图11 a 为实测数据反演网格结果，反演节点数为62 841 个，四面体单元数为 390 828 个，三角单元数为 783 079 个，边数为 455 091 个；图 11b为反演结果横切图，切片深度分别为10 m，图 8 设置磁化率模型Fig.8 Setting of the magnetic susceptibility modela俯视图；b立体图。图 9复杂模型反演结果Fig.9 Inversion results of complex modela异常体对应深度横切片图；bX和 Y方向纵切片图。839铀 矿 地 质第 39 卷-40 m，-70 m，-110 m 。如图 11所示，应用本文提出算法对实测数据进行处理得到较好的磁化率分布特征，模型边界分辨率高，异常特征明显，基于非结构网格其适用性强，具有模拟复杂结构的能力，对实测数据地形有较好的拟合。从实测数据反演 RMS误差曲线可以看出（图 12），随着迭代次数 RMS整体呈下降趋势，说明实测数据反演同样稳定收敛。基于以上试算说明应用本文提出算法进行大规模带地形磁测数据反演具有有效性，为大区域铀矿地质勘查提供一定技术支持。图 11 实测数据反演结果Fig.11 Inversion results of measured dataa实测数据反演结果；b反演结果深度500 m切片图；c反演结果深度1 500 m切片图；d反演结果深度2 000 m切片图。图 12 实测数据反演 RMS误差曲线Fig.12 RMS error curve of inversion result from measured data4结论本文开展了基于非结构化网格的大规模磁测数据三维正反演研究，分析了观测网格密度、边界节点数及局部网格细化等因素对正演精度的影响，探讨了利用灵敏度矩阵隐式存储的方法实现大规模磁测数据反演的算法，得到以下结论：1）磁异常的非结构化网格的有限单元正演在内存消耗和计算效率两方面均优于解析算法。2）非结构化网格有限单元正演的网格离散参数对计算精度有直接影响。提高观测网 840刘阳，等：基于非结构化网格的大规模磁测数据正则化反演研究第 5期格密度、增多边界节点数并采用局部网格细化技术，可提高磁测数据正演的精度和稳定性。3）非结构网格有限单元正演与灵敏度矩阵隐式存储的高斯-牛顿法反演相结合，有效减小了反演的内存消耗，降低了计算成本，适用于大规模三维磁测数据反演。参考文献1 程纪星，谢国发，乔宝强.音频大地电磁测深法与高精度磁法在相山铀矿田西部铀成矿有利远景预测中的应用 J.世界核地质科学，2013，30（2）：103-109.CHENGJixing，XIEGuofa，QIAOBaoqiang.Applicationofaudiofrequencymagnetotelluricsounding and high precision magnetic method in theprediction of uranium mineralization potential inwestern Xiangshan uranium field J.World NuclearGeology，2013，30（2）：103-109（in Chinese）.2 NABIGHIAN M N，GRAUCH V J S，HANSEN R O，etal.The historical development of the magnetic methodin exploration J.Geophysics，2005，70（6）：33-61.3 KAMINSKI V，HAMMACK R W，HARBERT W，et al.Geophysical helicopter-based magnetic methods forlocating wells J.Geophysics，2018，83（5）：269-279.4DARIJANI M，FARQUHARSON C G，LELIVRE PG.Joint and constrained inversion of magnetic andgravity data：A case history from the Mcarthur RiverArea，Canada J.Geophysics，2021，86（2）：79-95.5 DAI S，ZHAO D，WANG S，et al.Three-dimensionalnumerical modeling of gravity and magnetic anomaly ina mixed space-wavenumber domain J.Geophysics，2019，84（4）：41-54.6 NAGY D.The gravitational attraction of a rightrectangular prismJ.Geophysics，1966，31（2）：362-371.7 SINGH S K，SABINA F J.Magnetic anomaly due to avertical right circular cylinder with arbitrary polarizationJ.Geophysics，1978，43（1）：173-178.8 HOLSTEIN H.Gravimagnetic similarity in anomalyformulas for uniform polyhedra J.Geophysics，2002，67（4）：1126-1133.9 BARNETT C T.Theoretical modeling of the magneticand gravitational fields of an arbitrarily shaped three-dimensional body J.Geophysics，1976，41（6）：1353-1364.10FARQUHARSON C G，MOSHER C R W.Three-dimensional modelling of gravity data using finitedifferences J.Journal of Applied Geophysics，2009，68（3）：417-422.11LELIEVRE P G，OLDENBURG D W.Magnetic forwardmodelling and inversion for high susceptibilityJ.Geophysical Journal International，2006，166（1）：76-90.12JAHANDARI H，FARQUHARSON C G.Forwardmodeling of gravity data using finite-volume and finite-element methods on unstructured grids J.Geophysics，2013，78（3）：69-80.13OUYANG F，CHEN L.Iterative magnetic forwardmodeling for high susceptibility based on integralequationandGauss-fastFouriertransform J.Geophysics，2020，85（1）：1-13.14CAI Y，WANG C.Fast finite-element calculation ofgravity anomaly in complex geological regionsJ.Geophysical Journal International，2005，162（3）：696-708.15任政勇，汤井田.基于局部加密非结构化网格的三维电阻率法有限元数值模拟 J.地球物理学报，2009，52（10）：2627-2634.REN Zhengyong，TANG Jingtian.Finite elementnumerical simulation of 3D resistivity method based onlocally encrypted unstructured mesh J .Chinese Journalof Geophysics，2009，52（10）：2627-2634（in Chinese）.16ANSARI S，FARQUHARSON C G.3D finite-elementforward modeling of electromagnetic data using vectorand scalar potentials and unstructured gridsJ.Geophysics，2014，79（4）：149-165.17ZHOU F，TANG J T，REN Z Y，et al.A hybrid finite-element and integral-equation method for forwardmodeling of 3D controlled-source electromagneticinduction J.Applied Geophysics，2018，15（3-4）：536-544.18LI M，ZHANG Z，YANG J，et al.Integrated geophysicalstudy in the cemetery of Marquis of HaihunJ.Archaeological Prospection，2021，28（4）：453-465.19MAAG E，LI Y.Parameter selection workflow for adiscrete-valued gravity inversion with guided fuzzy c-means clustering C.Xi an：International Workshopand Gravity，Electrical&Magnetic Methods and TheirApplications，2019：444-447.20易柯，张志勇，李曼，等.基于 FCM聚类算法的二维RMT电阻率与介电常数联合反演 J.地球物理学报，2022，65（6）：2340-2350.841铀 矿 地 质第 39 卷YI Ke，ZHANG Zhiyong，LI Man，et al.Joint inversionof resistivity and permittivity of two-dimensional RMTbased on FCM clustering algorithm J.Chinese Journalof Geophysics，2022，65（6）：2340-2350（in Chinese）.21ZHANG J，WANG C Y，SHI Y，et al.Three-dimensional crustal structure in central Taiwan fromgravity inversion with a parallel genetic algorithm J.Geophysics，2004，69（4）：917-924.22MAAG E，CAPRIOTTI J，LI Y，et al.3-D gravityinversion using the finite element methodC/SEGTechnical Program Expanded Abstracts.Tulsa：Societyof Exploration Geophysicists，2017：1713-1717.23PILKINGTON M.3D magnetic data-space inversionwith sparseness constraints J.Geophysics，2009，74（1）：7-15.24AVDEEVD，AVDEEVAA.3Dmagnetotelluricinversionusingalimited-memoryquasi-Newtonoptimization J.Geophysics，2009，74（3）：45-57.25CAO H，WANG K，WANG T，et al.Three-dimensionalmagnetotelluric axial anisotropic forward modeling andinversion J.Journal of Applied Geophysics，2018，153：75-89.26JAHANDARIH，FARQUHARSONCG.3-Dminimum-structure inv