谱激电法是上世纪八十年代以来发展起来的一种新的电法勘探分支，可以分为时间域谱激电和频率域谱激电(也称复电阻率法).近年来，国内外学者对谱激电的正反演研究取得了不少成果：刘崧(1998)系统地介绍了谱激电的方法原理以及视复电阻率的理论计算与模拟；罗延钟等(1987)研究了同时考虑激电效应和电磁效应时复电阻率的数值模拟计算；徐凯军(2007)对2.5维复电阻率电磁场正反演进行了系统的阐述；李勇等(2011)研究了基于异常复电位的CR法数值模拟；范翠松等(2012)对2.5维复电阻率反演进行了应用试验；国外学者Hönig和Tezkan(2007)对时间域激电数据的一维和二维cole-cole参数反演进行了研究.

从上述研究成果来看，谱激电反演研究虽然取得了一些成果，但是还有许多问题有待解决，尤其是时间域谱激电.笔者认为，在不考虑电磁耦合的情况下，复电阻率的正演计算与直流类电法是类似的；而对于时间域谱激电，Yuval和Oldenburg(1997)提出的从激电数据中获取cole-cole参数的方法是比较可行的，这个思路将cole-cole参数的获取建立在时间域电阻率和极化率的反演的基础上，因此，要更好地研究谱激电正反演以及激电参数的反演，有必要对直流电阻率的正反演工作进行更为深入的研究.在2.5维直流电阻率正反演方面，徐世浙(1986)、阮百尧和徐世浙(1998)、万乐(2000)、汤井田等(2010)、潘克家等(2012)学者都做过深入的研究.近年来，任政勇和汤井田(2009)、汤井田和公劲喆(2010)、李长伟等(2010)学者对电阻率三维正演也做了研究.有限单元法是主要的模拟方法(熊彬和阮百尧，2002)，自1971年Coggon首先将有限单元法应用于直流电法数值模拟以来，其在地球物理中得到了广泛的研究和应用(周熙襄， 1983，罗延钟和孟永良，1986).

综上，2.5维直流电阻率有限单元正演的研究已经比较成熟，三维正演问题也在深入的研究中，本文在前人的研究基础上，对2.5维直流电阻率有限元正演问题进行了再次研究，主要是为了后续的时间域激电和复电阻率法的正反演研究服务，并对有限单元法求解偏微分方程有更深入的了解. 1 2.5维电阻率正演变分问题

在电阻率法中，点电源场是三维的，三维电场满足如下的Possion边值问题(潘克家等， 2012，万乐，2000)：

通常将具有一定走向的地质构造称为二维构造，其在走向上的电性特征可以认为是恒定的.假设二维地电断面的走向为y方向，则电导率在y方向上满足，因此二维构造的电性特征可以用一个二维坐标变量描述：σ=σ(x，z)，点电源场用一个三维变量描述：v=v(x，y，z).上述模型为二维、场源为三维的问题，称为2.5维问题.

2.5维问题的本质上仍是三维问题，使用有限单元法解决三维问题需要庞大的计算机内存(徐世浙，1994).上述2.5维问题，由于其对称性，可以借助傅立叶变换消除走向坐标y，将三维问题转换为带参数的二维问题(万乐，2000)，其波数域边值问题为

有限单元法求解偏微分方程，首先要将偏微分方程边值问题转化为对应的变分问题，与上式(2)对应的变分问题为

2 2.5维电阻率有限元正演

有限单元法的解题思路是对求解区离散化，将求解区剖分成若干个有限小单元，再对每个单元构造插值函数，用简化后的插值函数计算每个单元上的泛函，通过单元间的节点的函数值将各单元联系起来，获得整个求解区的泛函值，这样就将连续函数的泛函离散成了各单元节点上函数值的泛函，最后根据泛函求极值的条件，对泛函求一阶偏导数，得到各节点函数值满足的线性方程组，解线性方程组，得到各节点的函数值，从而得到有限单元法求解偏微分方程的结果. 2.1 有限单元求解

采用规则的矩形单元对求解区进行剖分，在求解区之外以等比步长向外扩展若干个节点来模拟无穷远边界，如图 1所示.

图 1

Fig. 1

图 1 网格剖分图 Fig. 1 The grid subdivision graph

为了对求解区进行更为精细的剖分和适应地形的变化，在每个矩形单元内再进行交叉对称三角形剖分(万乐，2000)，如图 2所示.

图 2

Fig. 2

图 2 网格单元示意图 Fig. 2 Sketch of finite element mesh

采用线性插值的方法对三角单元进行分析来求解变分问题，在三角形单元内，变分问题(3)为

将所有的三角单元的泛函相加，得到整个求解区的泛函

泛函求极值的条件是其对每一变量的一阶偏导数为零，即：

解上述线性方程组，便可以得到各节点的傅氏电位.为了节省内存，提高计算速度，本文将上述同一个矩形单元内的4个三角单元通过中心点的函数值建立联系，从而将各矩形中心点的函数值消去，降低了线性方程组的阶数(万乐，2000).

由于刚度矩阵是一个带状稀疏矩阵，对进行消元合并后的新的线性方程组的刚度矩阵采用了紧缩矩阵的存储方式，并使用不完全LU分解的稳定双共轭梯度算法求解线性方程组，得到线性方程组的 V .

2.2 反傅氏变换

通常，只需要对主剖面上的电位进行研究，由傅立叶逆变换得到

在实际计算时，由于只能利用有限个波数进行傅立叶逆变换的求解，因此得到的值是近似值，参与计算的波数越多，则求解的精度越高.如何选取适当的数值积分方法和波数，关系到傅立叶逆变换的求解精度，论文(万乐，2000)对此进行了详细的讨论，根据万乐博士的论文，通过下面的方法进行傅立叶逆变换：

总电位包含一次场电位和异常电位两部分，由于总电位在电源点附近是奇异的，因此，基于总场算法的有限单元法不能模拟点电源附近的总电位变化情况，而异常电位是不包含电源项的，因此基于异常电位的有限单元法能较好地模拟点电源附近的电位变化(徐世浙等， 1994，黄俊革，2003).

总电位v有两部分组成：

结合式(10)将式(11)分解得到异常场和一次场两部分：

异常场算法可以用下面的过程来实现：

求解同一点电源下的一次场傅氏电位和总场电位可以分别用如下线性方程组表示：

即：

设计如下两个理论模型，来考察异常场算法的优越性.

模型1：为一D型层状介质模型，如图 3所示，图 4为三极装置视电阻率正演结果.

图 3

Fig. 3

图 3 两层地电模型 Fig. 3 Schematic diagram of two-layer geoelectric section

图 4

Fig. 4

图 4 模型一正演结果 Fig. 4 Forward modeling results of modle 1

经过计算，总场算法的最大相对误差为37.71%，而异常场算法的最大相对误差仅为4.51%.

模型2： 如图 5所示二维低阻模型，图 6a为总场算法的二极装置视电阻率拟断面图，图 6b为异常场算法的视电阻率拟断面图.

图 5

Fig. 5

图 5 二维模型示意图 Fig. 5 Schematic diagram of 2D

图 6

Fig. 6

图 6 模型二正演结果 Fig. 6 Forward modeling results of modle 2

通过比较，使用总场算法时，在拟断面图的浅部有明显的低阻现象，而使用异常场算法时则不会出现这种情况.

通过以上两个理论模型的试验可以知道，由于靠近点电源处的奇异性，在点电源附近时，总场算法的正演结果是低于正确值的，模拟结果误差较大，甚至是错误，使用异常场算法则可以改善这种问题，它可以很好地模拟点电源附近电位变化情况. 4 模型验证

为了验证本文算法和程序的正确性，本文设计了如图 7所示的两个理论模型，使用本文的正演程序进行了正演模拟，并与国外某商业正演软件的模拟结果进行了对比.

图 7

Fig. 7

图 7 理论模型 Fig. 7 The theoretical model

对以上模型，采用偶极装置进行正演.图 8是采用本文的正演程序计算的结果，图 9是使用某软件正演的结果，通过对比，证明了本文所采用的正演方法和程序的编写是正确的.

图 8

Fig. 8

图 8 基于异常场的正演结果 Fig. 8 The Forward modeling results based on the anomaly electric field method

图 9

Fig. 9

图 9 某商业软件正演结果 Fig. 9 The Forward modeling results of a software

本文使用有限元法对2.5维电阻率进行了正演模拟，通过理论模型的验证，证明了本文的算法可靠，编写的程序是正确的，为后续进一步的正反演工作奠定了基础.

(1)在有限单元正演计算时，采用矩形网格和三角形网格相结合的网格剖分方法，在提高计算精度的同时，也减少了计算量，提高了效率.

(2)在正演过程中，采用了基于异常场的算法，提高了点电源附近的模拟精度.

(3)在解大型线性方程组时，对刚度矩阵采用了紧缩矩阵的存储方式，并采用不完全LU分解的稳定双共轭梯度算法，节省了计算机内存，提高了计算效率.