《文章编号1004-4140(二零一六)01-0047-08》.doc

文章编号：1004-4140（2009）01-0047-08 复杂三维大地电磁的有限单元法正演模拟( 张三四，李四五?，赵六七 （XX大学XX信息物理工程学院，长沙123456） 摘要：复杂二维和三维大地电磁模型的正演数值模拟具有一定的挑战性。对于复杂的二维和三维大地电磁正演问题，我们采用有限单元法进行求解。有限单元法最后形成一个线性方程组，系数矩阵是大型稀疏的带状对称复系数矩阵，并且其条件数远大于1，为严重病态矩阵，求解其对应方程组会遇到很多困难。不完全LU分解处理的Bi-CGSTAB迭代方法可用于该线性方程组的求解，并且具有速度快、精度高和稳定性好等优点；为了模拟无穷远边界及满足计算机的内存需求，在保证计算精度的情况下设计了非均匀网格剖分；在程序编制中，只存储有限元系数矩阵的非零元素，大大减少了正演计算的时间。通过对二维和三维模型电磁响应的计算，验证了算法的正确性。 关键词：大地电磁；有限单元法；正演模拟；Bi-CGSTAB算法 中图分类号：P631.3　　　文献标识码：A 二维/三维地电模型的大地电磁响应计算必须依靠数值模拟方法，数值解通常通过微分或积分方程的逼近和线性方程的求解来得到。对于大地电磁测深的正演模拟，最常用的3种数值模拟方法有：积分方程法、有限差分法和有限元法。总的来说，3种数值模拟方法在大地电磁测深正演计算方面都有一定的优势和不足：有限差分法的优点是方法简便易算，其缺点是，有限单元法与前述方法相比，在大地电磁测深正演模拟方面有独到的优势[1]。 1971年Coggon[2]首先将有限单元法应用在电磁法正演模拟中，他从电磁场能量最小原理出发，实现了二维地电断面有限单元法正演计算，不过由于有限单元法网格缺乏通用性，计算精度和速度未能达到实用水平；1976年William等[3]发展了有限单元法的剖分方法，采用矩形网格剖分，以解决二维大地电磁测深正演问题，使有限单元法向前发展了一步黄临平等[11]提出有限单元法的三维大地电磁模拟；王若等[12]用有限单元法正演模拟了三维CSAMT。 1　二维/三维大地电磁正演分析 由电磁场的Maxwell方程，角频率为（时间因子为）的电磁场的基本方程为： （1） 这里，E和H分别为电场强度和磁场强度，为电导率，为磁导率。 其中AB表示上边界，CD表示下边界。对于有限单元法求解三维大地电磁正演问题，由广义变分原理建立泛函得： （5） 有限单元法最后形成的总体刚度矩阵K是一大型稀疏、带状、对称复系数矩阵，并且其条件数远大于1，为严重病态矩阵，求解其对应方程组会遇到很多困难。对于二维大地电磁正演问题，K正定；而三维正演问题，K不正定。 2　正演计算策略 2.1　网格剖分分析 用有限单元法求解二维/三维大地电磁测深正演问题，首先要对研究区域进行剖分，节点与单元遍及整个研究区域。为了解决这一矛盾，我们将整个区域划分为2个区域：目标区域和网格外延区域，如图1和图2所示，目标区域为地质体赋存区域，也是数据的采集区域，以均匀网格剖分，网格外延区域的网格步长按大于1的倍数递增，在保证计算精度的情况下，减少网格剖分数，节省计算时间。 图2　三维模型不均匀网格剖分Fig.2　Non-uniform grid for 3D model 2.2　系数矩阵存储 用有限单元法做数值模拟时，由小单元集成所有节点的总体系数矩阵是一个大型稀疏矩阵，包含了大量的零元素，网格节点越多，零元素也越多。这里，我们采用matlab7.0为编程语言，对系数矩阵采取稀疏存储方式，即只存储系数矩阵的非零元素。 2.3　方程组求解 大地电磁二维/三维有限元正演计算最终归结为求解大型稀疏、对称，病态的复系数线性方程组。关于这类方程的求解国内外有较多的研究，直接解法有CS法、奇异值分解、LDLT分解和其改进方法等；迭代法有G-S法、Newton法、共轭梯度（CG）法及其改进方法。这里，我们采用不完全LU分解预处理的BICGSTAB算法[13]求解该线性方程组，该方法具有速度快、精度高和稳定性好等优点。 表3　人工神经网络训练与预测值Table 1　The analysis table of the sensitivity and specificity of DP diagnosis 试验编号 极间距 d∕mm 气压 p∕Pa 吸收率κ∕％ 渗层厚度δ∕μm 元素总质量分数w∕％ 试验值 预测值 试验值 预测值 试验值 预测值 1 1 1 70.900 70.587 34.5 34.579 87.496 87.437 2 2 2 61.200 60.871 36.5 36.380 89.