１引言

大地电磁法（ＭＴ）已广泛应用于地球物理探测各个领域．我国的山区面积比较大，在起伏地形条件下开展ＭＴ工作时，地形的影响产生电磁场畸变，这给资料处理和解释工作带来很大困难，因此有必要研究带地形的ＭＴ数值模拟算法．目前，在二维起伏地形和地电模型条件下，国内外学者采用过多种有限元剖分方式的二维ＭＴ有限元数值模拟．在有限元电性参数（如电导率、传播系数等）分块均匀的前提下，Ｗａｎｎａｍａｋｅｒ［１］、Ｍａｕｒｉｅｌｌｏ等［２］采用矩形三角网格剖分、双线性插值有限元进行数值模拟，陈小斌［３－４］、谢飞［５］用自适应地形三角网格剖分、双线性插值有限元数值模拟，周熙襄［６］、王绪本［７］、刘云等［８］采用自适应地形四边形网格剖分、双二次插值有限元数值模拟，对斜坡型地形产生的地形影响做模拟和分析；认为电性参数分块均匀的模拟方式，在单元块内的电性参数分布是均匀的，单元块与单元块之间的电性参数分布则不一定是连续的．但是在野外实际中，岩石、矿物的电性参数往往是变化的，因此采用电性参数分块连续变化的数值模拟方法，更符合野外实际．徐世浙、李予国等［９－１０］采用矩形网格剖分，进行电性参数分块连续变化、双线性插值、双二次插值等有限元数值模拟；对于矩形网格剖分，如果在有限元矩形单元块内或单元节点处填入空气介质，从而形成台阶型地形模拟地形影响［１１－１２］，是有其局限性的［４，１３］．阮百尧［１４］在起伏地形二维直流电阻率法正演数值模拟中，采用矩形内三角网格剖分、电导率分块连续变化、双线性插值有限元数值模拟．本文在此基础上，导出了起伏地形三角单元网格剖分电性参数分块连续变化的二维ＭＴ有限元数值模拟方法．

２二维ＭＴ变分问题

如图１所示，假设三维地质模型的电、磁等参数沿构造方向（ｘ方向）上无变化，这时可将三维模型问题转化为二维模型问题进行处理．此时，与二维ＭＴ边值问题相应的变分问题［１３］是其中，σ和μ分别为介质电导率和磁导率，ω为圆频率，ε为介电常数．下边界ＣＤ处离目标区域足够远，这时σ为ＣＤ处均匀介质电导率；左、右边界ＡＣ、ＢＤ离目标区域足够远．３有限单元法用有限单元法求解式（１）的二维ＭＴ变分问题，具体步骤如下：

３．１网格剖分对于ＴＥ模式，在高频时，空气介质中位移电流的影响不可忽略［３］，因此研究区域为空气和地下；对于ＴＭ模式，场值Ｈｘ在空气中近似为常数（在ｙ和ｚ方向的偏导数近似为零），即不考虑空气介质中位移电流的影响，此时ＡＢ边界设在地面上，研究区域为地下．如图１所示，研究区域Ω用三角单元进行网格剖分．如图１ａ所示，当模拟水平地形时，考虑到地下介质体形状的任意性，首先对研究区域进行矩形网格剖分，再在每一个矩形网格内剖分出４个三角单元网格．这样，一方面，避免了三角网格过于尖锐的情况；另一方面，可利用三角形的斜边模拟地形线和任意介质体倾斜的界面；如图１ｂ所示，当模拟起伏地形时，在图１ａ网格剖分的基础上，地空边界采用沿实际地形线进行网格剖分．

３．２线性插值如图２所示，假设三角单元ｅ内电磁场ｕ和电性参数τ，λ，ｋ呈线性变化，则在每一个单元中：

３．３单元分析在整个研究区域Ω中，将泛函（１）式离散化，表示为所有三角单元ｅ的线性组合，即

３．３．１单元分析１ｕ对ｙ求偏导数，有

３．３．２单元分
线积分∫ｌＮａ１Ｎｂ２ｄΓ＝ａ！ｂ！（ａ＋ｂ＋１）！ｌ，其中ａ，ｂ为非负整数，ｌ为单元底边的长度［１３］．对于（２）式中第三个积分项只对ＣＤ边界进行．当单元的２３边落在ＣＤ边界上时，线积分

３．４系数矩阵总体合成及求变分将研究区域Ω内每个单元系数矩阵ｋｅ按照总体节点号进行扩展，得到ｋｅ，相加得到总体系数矩阵Ｋ，即如图３所示，采用矩形内三角网格剖分方式，比较矩形网格剖分来说，增加了一个中间节点５，大大增加了计算量．但是节点５只与该矩形内４个角点有关系，而与其它不相邻的矩形网格没有任何直接联系，所以在单元分析中，运用高斯消元法事先可消去这个节点．节点５只是一个虚设的节点，并不包含在总体节点数中．到解方程结束后，亦可根据矩形网格内４个角点的已知场值和节点系数之间的关系，直接计算出节点５的场值．这样，一方面不增加节点总数，节省了计算量，另一方面每个三角单元仍具有

４辅助场、视电阻率和阻抗相位的计算

４．１辅助场的计算辅助场是通过计算主场值沿某方向的方向导数得到的［１７－１９］．当各节点主场值求出后，由３．２节可知，单元内场值ｕ可表示为ｕ＝∑３ｉ＝１Ｎｉｕｉ，则ｕ对ｙ求偏导数，得元边上任意点或三角单元的节点为实际测量点，此时要考虑到相邻单元间场值偏导数的不连续性问题［１３］．对于彼此相邻的单元、相同单元边或节点上场值的偏导数，可以在所属不同单元中予以分别计算，再取平均值．本文以地面三角单元的节点为实际测量点．不难发现，对于ＴＥ模式，需要计算同一测点的地表和地表上一空气层所有单元的场值偏导数，再取平均值；而对于ＴＭ模式，则只需求出同一测点的地表所有单元的场值偏导数，再取平均值即可．

４．２视电阻率定义和阻抗相位的计算对于ＴＥ模式：ｕ＝Ｅｘ，辅助场Ｈｙ＝τ?Ｅｘ?ｚ．在二维起伏地形条件下，Ｅｘ、Ｈｙ在野外实际中均可测５模型计算为验证本文正演算法的有效性，首先与水平层状一维连续介质理论模型的解析法、以及前人已有地形模型做计算分析，之后用该算法对倾斜界面异常体模型做计算分析．在以下所有模型的计算中，有限元网格的剖分方法如下：如图１所示，在横向网格中，目标区域采用等间隔网格，间隔距等于点距．左、右稀疏网格数各为１８个，网格大小等于点距乘以系数，其中１８个系数依次分别为１，１，１，１，１，２，２，２，２，４，４，４，８，８，１６，３２，６４，１２８；在纵向网格中，空气稀疏网格坐标（ＴＭ模式无空气网格，单位为ｍ）１４个，依次分别为１０００００，５００００，１００００，５０００，１０００，５００，２００，１００，５５，３５，２５，１５，１０，５．地下网格的剖分，从１～１０层的网格大小（单位为ｍ）依次分别为５，５，５，５，１０，１０，１０，１５，１５，２０，第１层网格大小的默认值为５ｍ．当然，也可以任意设置第１层网格大小，这样，从１～１０层的网格大小则相应地等比例放大或缩小．从１１层及以下的网格剖分为等比网格，该层的网格大小等于上一层网格大小乘以等比系数，等比系数一般为１．１或１．２．这样，横向网格数则取决于测点数的多少，纵向网格数则取决于探测深度的选择；当模拟起伏地形时，上、下地表附近还要加上地形网格，地形网格数等于所有测点的高程数（对于重复的高程，只计１个），地形网格大小等于相邻高程之差．

５．１一维连续介质模型如图４所示，三层一维连续介质地电断面模型，第一层电阻率为１００Ωｍ、层厚度为１０００ｍ；第二层为连续介质层，介质电阻率随深度线性减小，其变化关系如图中所示，层厚度为１０００ｍ；第三层电阻率为１Ωｍ；测点区域为０～４ｋｍ，点距为１００ｍ，共４１个测点，频率范围为１０３～１０－３Ｈｚ，以１０为底的对数间隔采样，共５１个频点．在一维层状介质的解析法中，将中间的连续介质层剖分为１０层进行模拟．在二维有限元数值模拟中，取测线中心处第２１号点测点为研究对象．则解析法与二维数值解法的结果对比如图５所示．在ＴＥ、ＴＭ视电阻率曲线图和阻抗相位曲线图上，解析法和数值解曲线形态、值的大小基本一致，这表明本文方法能够有效地模拟连续介质模型．正演结果的数据统计表明，视电阻率、阻抗相位的解析法结果和数值模拟结果均方根误差均小于１％．本文的数值模拟结果与解析法结果之间的误差分析主要为：一方面，与有限元数值模拟的网格剖分有关系，如高频网格和低频网格的剖分方式是不一样的；另一方面，解析法是基于电性参数分层均匀的计算方式，而本文的数值模拟方法是基于电性参数分块连续变化的模拟．

５．２地形模型如图６下部所示为圆柱凹陷地形模型，Ｗａｎｎａｍａｋｅｒ等［１］和徐世浙等［２０］在文献中对这个模型做过模拟．测点区域为０～０．４ｋｍ，点距为１０ｍ，共４１个测点，圆半径为５０ｍ，频率为０．０１Ｈｚ，计算ＴＭ模式下视电阻率．本文数值模拟结果如图６上部曲线图所示，对照文献［１，２０］中模拟结果，在０．０１Ｈｚ频点处，ＴＭ模式下，三种模拟方法与地形测点对应的视电阻率起伏形态和值的大小基本一致．图７所示为山脊地形模型，在文献［１，２０］中对这个模型做过模拟．测点区域为０～４ｋｍ，点距为５０ｍ，共８１个测点，起伏落差为１００ｍ，山底宽度为２．４ｋｍ，频率为１０Ｈｚ．本文数值模拟结果如图８所示．对照文献［１，２０］中的模拟结果，在１０Ｈｚ频点处，ＴＥ、ＴＭ模式下，三种模拟方法与地形测点对应的视电阻率、阻抗相位起伏形态和值的大小基本一致．通过对以上两例均匀大地起伏地形模型的模拟，得到了与前人研究相符的结论．分块均匀法认为地面或地下单元内节点之间的电性参数是均匀分布的，而本文的分块连续法则认为地面或地下单元内节点之间的电性参数是连续变化的（空气介质除外）．显然，分块均匀法是分块连续法的特殊形表现形式，分块连续法更具有普遍性．用三角网格模拟起伏地形、分块连续变化的剖分方式且更加符合野外实际地质情况．

５．３倾斜界面异常体模型如图９所示为“Ｗ”字形倾斜界面异常体模型，背景电阻率为１００Ωｍ，异常体电阻率为２０Ωｍ，其上顶面距地面１００ｍ，下底面距地面５００ｍ，异常体上顶面外宽４８００ｍ，下底面外宽２４００ｍ．测点区域为０～８ｋｍ，点距为１００ｍ，共８１个测点，频率为１０３～１０－１Ｈｚ，以１０为底的对数间隔采样，共４１个频点．有限元数值模拟网格剖分为：最大探测深度选择为１０ｋｍ，等比网格系数选择为１．１；ＴＥ模式，网格数为１１６×６６（横向网格数为１１６，左、右稀疏网格数各为１８；纵向网格数为６６，其中包括空气稀疏网格数为１４）；ＴＭ模式，网格数为１１６×５２（横向网格数同ＴＥ模式；纵向网格数为５２，没有空气网格）．大型稀疏矩阵的存储方式均按文献［１３，１５］中变带宽方式存储，用不带平方根的Ｃｈｏｌｅｓｋｙ分解法求解线性方程组［１３，１６］．本文数值模拟方法的结果剖面图如图１０所示，计算每一个频点的时间为０．７８ｓ（ＴＥ和ＴＭ模式同时计算）．在ＴＥ、ＴＭ两种极化模式下，视电阻率、阻抗相位剖面上都出现明显低阻异常，“Ｗ”形异常形态清晰可见，表明采用三角网格模拟任意倾斜界面异常体的方法是有效的．

６结论

本文提出三角单元网格剖分、电性参数分块连续变化二维ＭＴ有限元数值模拟方法，比较矩形单元剖分、电性参数分块均匀的模拟方法，这样既不增加节点，节省了计算时间，又能很好模拟野外实际地形和任意地电体模型．通过对一维连续介质模型模拟计算，与解析法计算结果对比，验证了本文算法的可靠性；对前人地形模型模拟计算，结果与前人结果相符；对任意倾斜界面异常体模型进行模拟，结果表明本方法是有效的．同时指出，本方法仍然是考虑各向同性介质情况，对各向异性介质分块连续变化的数值模拟方法将是本文的后续研究工作．

编辑提示，此文是论文格式网为朋友们总结并提醒职称及职称考试的相关事项。希望对朋友们有所帮助。