张胜业，孟永良，黄理善，石砚斌，王长春，曾歌明

张胜业，孟永良，黄理善，石砚斌，王长春，曾歌明_大人地球物理科学

张胜业¹，孟永良¹，黄理善²，石砚斌³，王长春⁴，曾歌明⁵

（1.中国地质大学地球物理与空间信息学院，武汉，430074；2.桂林矿产地质研究所，桂林，541004；

3.防灾科技学院，三河市，065201；4.山西省煤炭工业厅煤炭资源地质局，太原，030045；

5.中船重工集团第722研究所，武汉，430070）

摘　要：以稳定电流场中2.5维地电模型电阻率计算为基础，应用Cole-Cole模型，把极化效应的影响考虑到电阻率的计算中，应用有限单元法及边界单元法对2.5维起伏地形地电模型的电阻率和极化率进行频谱和时谱激电资料的正反演方法研究，为谱激电资料的二维反演提供一个实用的方法。在实验室使用GDP32II多功能电法仪和国产DJS-8直流激电仪等仪器对人工制作的几十块标本和在矿区采集的天然岩矿石标本进行了复电阻率的频谱和时谱测量，获得人工合成岩石、矿石标本和天然标本的复电阻率的频谱和时谱曲线、极化率和电阻率等参数。采用较成熟的Cole-Cole理论模型拟合标本的复电阻率频谱和时谱曲线，获得标本的4个Cole-Cole模型参数。通过研究岩石、矿石的4个Cole-Cole模型参数随其导电矿物的成分、结构、含量和颗粒大小等因素的变化规律，从中总结出一些有意义的研究成果。

关键词：谱激电法，2.5维，有限单元法，边界元法，模型实验Cole-Cole模型

前　言

目前谱激电法被广泛用于环境监测、矿产勘探和水文地质学领域，具有获取地电信息丰富、抗干扰能力强等优点，是目前常用的电法勘探方法之一^［1～12］。由于三维电法模拟计算量大，以目前计算机的计算速度来看，实用价值不高，因此，2.5维问题被广泛地研究^{［13～19］}，2.5维问题即三维场源，二维地质体。通过在地质体走向方向上进行傅立叶变换可以将三维问题转变为波数域的二维问题，其计算量远远小于纯三维问题。笔者采用有限单元法，研究在地形不平条件下，不考虑电磁干扰的频谱和时谱激电2.5维的正演模拟和反演算法^{［20～22］}。

岩石、矿石的复电阻率振幅谱和相位谱包含有岩石、矿石导电矿物的成分、结构、含量和颗粒大小等丰富的信息。系统地研究人工合成标本和天然标本的岩石、矿石的复电阻率频谱和时谱特性，对详细评价异常体具有重要的理论和实用价值^{［23～26］}。

实验用水泥、河砂、石墨、黄铜和铁块等材料人工合成导电矿物成分、结构、含量和颗粒大小不同的多相岩石、矿石标本和在矿区采集的天然标本，使用GDP32II多功能电法仪的岩心测试设备、国产DJS-8直流激电仪、电压表和电流表等仪器和仪表，获得人工合成岩石、矿石标本和天然标本的复电阻率的频谱和时谱曲线、极化率和电阻率等参数。采用较成熟的Cole-Cole理论模型拟合标本的复电阻率频谱和时谱曲线，获得标本的4个Cole-Cole模型参数。研究岩石、矿石的4个Cole-Cole模型参数随其导电矿物的成分、结构、含量和颗粒大小等因素变化而变化的规律，从中总结出一些有意义的研究成果。

1　2.5维起伏地形模型的频谱激电法的正反演研究

1.1　2.5维地形起伏模型的频谱激电的正演研究

1.1.1　方法原理

以稳定电流场中2.5维地电模型电阻率计算为基础，应用Cole-Cole模型，把极化效应的影响考虑到电阻率的计算中，应用有限单元法及边界单元法对2.5维起伏地形地电模型的电阻率和极化率进行正演计算^［20］。

1.1.2　有限单元法的变分方程

对于点源电流地下电场满足的方程：

为了消除沿走向的坐标y（使变为二维问题），对上式两端分别作余弦傅氏变换，变换后变为：

式中：σ＝σ（x，z），V＝V（λ，x，z），λ为傅氏变换变量或称为波数。其中V为电位U的傅氏变换。

以上在波数域的电位方程等价的变分问题为：

式中：

1.1.3　网格剖分

应用如图1所示的网格剖分方法，在以往矩形剖分的基础上加入地形的影响成为平行四边形剖分，这样计算刚度矩阵的公式就发生了变化。但是，求得的刚度矩阵仍然是一个正定的大型稀疏带状矩阵，可以应用Cholesky分解法求解方程组。

1.1.4　边界元法应用

在计算转换电位时应用二次场算法，既不显著增加计算量，又大大提高了模拟的精度。在地形不平的情况下一次场的计算没有解析解，因此应用精度较高的边界单元法来求解此时的一次场。

边界单元法只在研究区域的边界上剖分单元，这样可以降低求解问题的维数，使最后形成的代数方程组的元数减少，经常用于二维直流电阻率法的地形改正。借助格林公式及二维介质霍姆赫兹方程的基本解，即可把转换电位V满足的霍姆赫兹方程及边界条件等价地归化为如下的边界积分方程：

图1　网格剖分示意图

式中：C_i是常数，为第二类修正一阶贝塞尔函数。

把地形边界Γ分割成N-1个单元，N个节点。在每一个单元ΔΓ_j上取插值函数V为一常数，在ΔΓ_j上每一点的值都相等，此时式（4）可写形式为：式中：积分∫Δ_Γjq^＊（x，y）dΓ利用5点高斯—勒让德求积公式求解，上式即为关于V的线性方程组，可求得V的值。

1.1.5　程序验证

为了验证算法的正确性，使用一个一维地电模型视电阻率的解析解进行对比。所用模型为一个水平层状介质模型，其参数分别为ρ₁＝100Ω·m，h₁＝10m，ρ₂＝500Ω·m，h₂＝10m，ρ3＝50Ω·m。对该地电模型使用AMN三极装置进行测量，解析计算和有限元正演模拟结果的对比如图2所示。

由计算结果可以看出，本文所用算法精确度较高，有限单元法结果和解析计算结果的相对均方误差仅有0.73%。

1.2　2.5维带地形模型的频谱激电的反演研究

1.2.1　反演方法

在地表进行频谱测量工作时，每个频率每次供电，都可以得到一组数据，则Nd次供电测量得到Nd组数据，分别为ρ_s^R₁，ρ_s^R₂，…，ρ_s^R_n，它们一起构成原始测量数据｛ρ_s^R｝。对于一个地电断面模型，我们给定一组初始模型参数σ^o₁，σ^o₂，…，σ^o_m，它们构成初始模型向量｛σ^o｝。使用这些模型参数，我们通过Nd次正演计算（不同供电点位置），可以获得Nd组理论计算结果，分别为ρ_s^C₁，ρ_s^C₂，…，ρ_s^C_n，它们一起构成理论计算数据｛ρ_s^C｝，即｛ρ_s^C｝是｛σ^o｝的函数。我们要对高密度电法的测量数据进行反演，就是要获得地电断面的电导率分布，即求出σ₁，σ₂，…，σ_m，使得我们模型计算结果与野外测量数据结果吻合最好。从数学上来说，就是要使：

图2　解析计算和有限单元法模拟结果对比图

i＝1

i＝1

要解式（6），通常要解下列法方程：

式中：Δρ_s＝｛Δρ_si｝＝｛ρ_s^R_i-ρ_s^C_i（σ^o）｝（i＝1，2，…，n），ρ_s^R_i为野外实测数据，ρ_s^C_i（σ^o）为使用模型参数｛σ^o｝作正演计算获得的结果；A称为雅可比（Jacobi）矩阵，其元素为：

Δσ称为模型修正量，Δσ＝｛Δσ_j｝＝｛σ_j-σ^o_j｝，（j＝1，2，…，m），｛σ^o_j｝是我们给定的初始模型参数，而｛σ_j｝就正是我们要求的模型参数。即：

式（7）是一个高阶线性方程组，获得雅可比矩阵A后，可以通过矩阵分解求逆来解此方程，以求得Δσ：

或者直接求A的广义逆矩阵A，来求解此方程：

此线性方程组的解法已比较成熟，应用也很广泛，在此，我们不详细讨论式（10）或式（11）的解法。

1.2.2　反演算例

图3为起伏地形下水平低阻体的模型示意图，模型参数如图中所示。对下图模型采用类似高密度电法方式，对二极装置进行计算，其中，电极个数为40，电极距为5m，最小排列系数为3，最大排列系数为18，供电频率为2⁰～2³ Hz 4个频段，根据前面介绍的反演方法进行反演计算。

图3　起伏地形下低阻体模型图

ρ_s0和m_s分别经过5次迭代计算，得到计算结果如图4和图5所示，数据拟合均方相对误差分别为4.25%和9.62%。在对起伏地形模型反演时可以看出，地形因素仅对视电阻率的影响较大，对其他激电参数影响较小。

图4　反演得到的频率为0时的电阻率ρ0模型图

图5　反演得到充电率m模型图

由以上算例可以看出，反演的结果已经可以满足一般地质要求。对于起伏地形的反演，基本消除了地形影响，并很好地反映出了位置、形状以及电性。

1.3　小结

（1）频谱激电法采用类似高密度电法装置测量，可以保证获得广泛的地下电性参数信息。

（2）在进行正演计算时，采用二次场算法，由于其计算的是由地下不均匀体产生的二次场变化，因此计算精确度要高于总场算法。

（3）在进行视参数到真参数的反演时，加入光滑度矩阵，可以使得反演结果较为光滑并接近于预设的基本模型。反演充电率m_s时，直接使用反演视电阻率ρ_s的雅克比矩阵，不需要额外的计算，从而提高了运算速度。

（5）在对测量记录的视复电阻率谱反演Cole-Cole模型视参数时，根据Cole-Cole模型公式，采用递推的方法进行反演，方法简单有效，计算速度快，不需要进行参数估计，没有额外的变量，算法稳定。

2　2.5维起伏地形模型的时间谱激电法的正反演研究^{［21，22］}

2.1　2.5维地形起伏模型的频谱激电的正演研究

由于其方法原理和计算方法与频谱激电法类似，所以，这里就略去了。

2.1.1　时域电阻率的Cole-Cole表达式

通过对频率域的Cole-Cole模型进行拉氏变换，可以得到Cole-Cole模型的时间域表达式：

上式可以通过Guptasarma数值滤波算法［26］近似求得。

然后利用“等效电阻率”，应用各种电阻率数值模拟方法求得地表的视电阻率，进而求得某一断电时间t下的视极化率：

2.1.2　程序验证

对于验证所编程序的正确性，需要应用有解析解的模型，以及对前人所做的工作进行验证，下面将从水平层状大地模型对所编的程序进行验证。假设地下介质模型为水平3层的H型模型，第一层厚度为5m，电阻率为1 000Ω·m，第二层厚度为5m，电阻率为100Ω·m，第三层电阻率为1 000Ω·m，厚度无限大。得到的模拟结果如下。

在图6中，红色曲线为数值模拟结果，绿色曲线为没有应用二次场算法的数值模拟解，蓝色曲线为应用了二次场算法后的数值模拟解，其中总场算法平均相对误差为4.13%；二次场算法平均相对误差为3.03%。从上可以看出，在本文中所应用的有限单元模拟的程序是正确的，并且精度很高，尤其是应用了二次场算法后，大大提高了较小极距时数值模拟的精度。

2.1.3　正演算例

图7所示的是山谷地形下直立板状极化模型正演模拟结果图。

2.2　2.5维带地形模型的时间谱激电的反演研究

笔者主要采用类似于高密度电法的装置对所得到的视电阻率值进行反演，求得每个时间地下介质单元的“真”电阻率值。在实际的反演中，电阻率的反演方法有很多种，在文中笔者采用了应用比较广泛的广义逆法。在电阻率的多种反演方法中，雅克比矩阵的计算都是不可避免的，也是电阻率反演中的一个最重要的步骤已经过国内外学者的研究，Tripp等^［27］提出了应用电位函数与电导率参数的关系，求取雅克比矩阵的方法，大大降低了求取偏导数矩阵时的计算量。国内，阮百尧^{［28，29］}、万乐^［30］、黄俊革^［31］等也在这方面做了大量的工作，并且应用互换定律，只需做一次正演计算，就可以得到雅克比矩阵的值，进一步提高了计算的速度。

图6　三层H型模型模拟结果对比图

图7　山谷地形直立板状极化模型正演模拟结果示意图

（a）模型示意图；（b）中梯装置视极化率曲线；（c）二级装置视电阻率拟断面图；（d）二级装置视极化率拟断面图

由Seigel^［32］的研究表明，视极化率和地下网格的“真”极化率的关系可以表示成

m_s＝-A·m

式中：A即为上面提到的视电阻率对单元电导率的偏导数矩阵，因此在求解得到雅克比矩阵后可以容易的反演得到地下单元的极化率值。

所设的反演模型如图8所示，为一山脊地形下低阻高极化体，模型背景（图中用蓝色表示部分）电导率为0.01S/m，极化率为0.01，频率相关系数为0.25，时间常数为1s；异常体（图中用红色表示部分）电导率为0.1S/m，极化率为0.5，频率相关系数为0.45，时间常数为3s。

图8　山脊地形下低阻高极化体模型示意图

电导率反演结果如图9所示。

图9　电导率模型反演结果

极化率反演结果如图10所示。

图10　极化率模型反演结果

2.3　小结

（1）本文中应用平行四边形网格对带地形模型进行剖分，在前人矩形剖分的基础上重新计算了刚度矩阵，实现了有限元正演工作。

（2）在有限元正演工作中应用二次场算法，尤其在极距比较小的情况下，大大提高了模拟精度。

（3）反演过程中应用了电位函数与电导率参数的简单关系求取雅克比矩阵的方法，大大降低了求取偏导数矩阵时的计算量。反演中加入了固定的阻尼因子，抑制了雅克比矩阵中的小奇异值，降低了求取线性方程组时的病态效应，使得反演的稳定性有所提高。

（4）对理论模型进行反演计算，反演结果与已知模型基本吻合，达到了重建地下电性参数的目的。在反演极化率时应用了偏导数矩阵与视极化率的关系，直接对极化率进行反演。

3　多相矿化岩石标本模型复电阻率的物理模拟及研究^{［23，24］}

通过采集和制作多相介质岩石、矿石标本模型，在频率域和时间域研究多相介质模型结构参数变化和复电阻率参数的变化规律。

这项研究工作主要是2007年和2009年进行的，由本项目两位成员曾歌明和黄理善研究生完成，并分别将此项研究作为他们硕士论文的主要内容。多相矿化岩石标本的获取一般采用两种形式。

一是从矿区采集矿石标本，近两年分别到内蒙古白音诺尔铅锌矿，内蒙古金厂沟梁金矿和山西垣曲县胡家峪铜矿等矿区采集到铅锌矿、金矿和铜矿的矿石及围岩的标本。

二是多相矿化岩石标本模型的制作，采用水泥、石英砂、石墨粉（颗粒）、金属球颗粒和块状金属等来模拟围岩和矿物。标本为10cm×10cm×10cm的正方体，制作时不施加额外压力，共制作了40块人工合成标本，它们的编号分别为RG01，RG02，…，RG40。

3.1　人工标本的制作和天然标本的处理

采用碳酸盐水泥和干净的细河砂（石英砂）模拟围岩（脉石）物质，导电矿物有不锈钢钢珠、铁粉、铁块、黄铜（粉末状、颗粒状和条块状）、石墨等，导电矿物的颗粒大小、形状、连通程度等按不同要求改变。使用人工合成岩石、矿石来进行研究，制作花费时间短，材料费用低。

天然标本采自内蒙古赤峰市巴林左旗的白音诺尔铅锌矿矿区，岩石、矿石类型有大理岩、矽卡岩、闪长岩、次生石英岩、含闪锌矿矽卡岩和闪锌矿矿石等，所有标本均切割成边长为5cm的正方体，测试前用普通自来水浸泡48小时。

同时对采集的天然岩石、矿石标本进行切割加工处理，使其成为规则形状，方便实验测试。共制作了40块人工合成标本，它们的编号分别为RG01，RG02，…，RG40。

3.2　岩矿石标本模型时间域谱激电的物理模拟及研究^［23］

我们对采集的矿区岩矿石标本模型进行了时间域谱激电的物理模拟工作，实验所用的标本是在内蒙古赤峰市巴林左旗白音诺尔铅锌矿区所采的岩石和矿石标本。

3.2.1　实验设备及技术参数

实验设备及技术参数如下：

（1）实验中采用的标本采自内蒙古赤峰市白音诺尔铅锌矿矿区，包括不同岩性。

（2）实验场所为室内，实验在一标本盒内进行，盒长15cm、宽5cm、深5cm。

（3）采用装置为对称四级装置，供电电极距AB＝15cm，测量电极距MN＝5cm。

（4）发送仪器为NANOTEM NT-20发送仪。

（5）接收仪器为GDP-32II型多功能电法仪。

（6）工作电流固定在0.03A，工作电压为0～24V。

实验系统是利用GDP-32II中时间域激发极化方法中的直流激电测量模块（TDIP），双向短脉冲直流供电，用对称四极电极装置，在标本盒中对切割成正方体的岩石标本进行测量。试验中采用的标本采自内蒙古赤峰市白音诺尔铅锌矿矿区，经切割处理为大小5cm×5cm×5cm左右的正方体，样品岩性包括铅锌矿矿石、变质岩、岩浆岩等围岩，为保证结果的可靠性，所有测试标本均取自岩石露头上或矿坑中的典型样品，测量前在水中浸泡2天以上。

3.2.2　实验结果统计和分析

从本次实验结果来看，不同的岩矿类型的激电特征有各自的特点，观察测得的充电率（表1），可以发现不同类型的岩石标本衰减曲线差异很大，基本可以把矿石与围岩区分开来，但相同类型也存在差异。

表1　不同岩石类型充电率值的范围

一般围岩的充电率都远小于含矿标本，但是部分含矿的围岩（如表1中的矽卡岩和闪长岩）也会出现较大的值。

把实验得到的结果以岩石类型分类，选择了12个时窗测量值，测量时间从断电后，每类选择一组数据画出二次场充电率积分衰减曲线图，如图11所示。

从图11看到，闪锌矿的充电率最高，含方铅矿的闪锌矿次之，部分矽卡岩因含金属矿物充电率也比较高，其他一些围岩充电率则很低。

由图11可以看出，对第二个时窗（断电后0.453 125s开始），不同类型部分岩样标本的充电率时窗积分衰减曲线存在差异。充电率值最大的是闪锌矿矿石标本，所测值的范围也最宽，从90直到10左右。含方铅矿的矿石其次，值的范围从60直到10以下。矽卡岩中含矿的标本充电率比较大，可以达到15左右，而不含矿的值只有2.5左右。其他围岩，包括角岩、闪长岩、大理岩、次生石英岩、火山凝灰岩等都在2.5左右，没有明显区别。含矿标本最后一个时窗的值往往也较高，最高可以达到50以上，也有部分下降很快，接近于零，而非矿标本最后一个时窗充电率基本都接近于零。

图11　不同类型部分岩样标本的充电率时窗积分衰减曲线

3.2.3　实验数据的反演与分析

本次实验数据的反演结果是人工交互拟合的，基本步骤是先代入估计的参数初值，得到理论曲线、实测曲线和误差值，然后按得到的二次场衰减曲线与参数变化规律，用人机交互的方法进行拟合，当拟合结果的误差较小时反演结束，下面给出方铅矿的拟合结果图，如图12所示。

图12　方铅矿标本拟合结果示意图

从拟合结果（表2）可以看到，不同岩石类型标本的m、c、τ值的特征与前人所述基本一致，矿石的激电谱拟合得较好，矿石的m和τ值都相对较大。

表2　典型岩矿标本粗略拟合值

3.3　岩矿石标本模型频率域复电阻率的物理模拟及研究^{［24，25］}

实验用水泥、河砂、石墨、黄铜和铁块等材料人工合成导电矿物成分、结构、含量和颗粒大小不同的多相岩石、矿石标本，使用GDP32II多功能电法仪的岩心测试设备、国产DJS-8直流激电仪、电压表和电流表等仪器和仪表，获得人工合成岩石、矿石标本和天然标本的复电阻率振幅频谱、相位频谱、相对极化率和直流电阻率。采用Cole-Cole理论模型拟合标本的复电阻率振幅频谱曲线，反演获得标本的4个Cole-Cole模型参数。研究岩石、矿石的4个Cole-Cole模型参数随其导电矿物的成分、结构、含量和颗粒大小等因素变化而变化的规律。

3.3.1　人工标本与测量装置

人工合成标本模型如图13所示。标本复电阻率测量采用四极装置（如图14），当电流I（ω）通过标本时，在标本中部两个等位面之间的电位差ΔU（ω）与标本的电阻率ρ（ω）成正比：

式中：ω为角频率；K为装置系数，它与标本的尺寸，形状以及电极间的距离有关；ΔU（ω）为电位差（mV）；I（ω）为电流（mA）；ρ（ω）为复电阻率振幅（Ω·m）。

图13　人工合成标本模型简图

图14　实验用标本复电阻率测量装置

由于人工合成标本的尺寸为正六面体，并且AB的距离比标本的边长大，所以此处：

式中：S为标本的横截面积；L为测量电极M、N之间的距离。

实验使用的是ZONGE公司生产的GDP32-Ⅱ多功能电法仪及其相应的实验室岩心测试设备。为保证电流只能且均匀地通过标本，在测试槽与标本两侧接触处用绝缘橡皮泥隔开，槽内用自来水与标本的两端接触面保持良好接触。使用四极装置测量，发射A、B极用铜板电极，接受M、N极用Cu-CuSO₄不极化电极。

3.3.2　参数反演

前人的研究已经证实了大多数岩石、矿石的激电谱可以用Cole-Cole模型描述（Pelton，1978），Cole-Cole模型表达式为：

式中：ρ₀为零频率时的电阻率；m为极限极化率（充电率）；τ和c分别为时间常数和频率相关系数，它们是用来评价激电异常的参数。

根据式（16）用人机交互方法对复电阻率振幅曲线进行拟合，求得标本的4个Cole-Cole模型参数：零频时的电阻率ρ₀、极限极化率m、时间常数τ和频率相关系数c。

3.3.4　岩矿石复电阻率特性分析

1.导电矿物体积含量与标本Cole-Cole模型参数的关系

（1）浸染状标本。实验用黄铜粉、铁粉和石墨粉模拟了浸染状岩石、矿石［图13（a）］，铁粉和黄铜粉的颗粒尺寸小于0.2mm，石墨粉的颗粒尺寸小于0.05mm，图15为典型浸染状人工标本复电阻率测量结果，同类岩石、矿石按导电矿物含量不同进行分组对比。

水泥∶河砂＝1∶2；含黄铜粉8%参数反演结果：ρ₀＝10.1Ω·m，m＝0.36，τ＝3.5×10^-5s，c＝0.49

图15　RG22标本测量结果

测量结果表明，水泥和河砂含量均为50%的RG1标本，不含良导电矿物成分，其极限极化率m与时间常数τ均较低；在测试频率范围内，其相位绝对值在0～30mrad之间，极化效应较弱。含粉粒状铁粉的人工标本RG5、RG4和RG3，铁粉含量分别为5%、10%和20%，其极限极化率m值随铁粉含量增大而明显增大；时间常数τ也随铁粉含量增大而增大，但其τ值均较小，只有10^-5s数量级；频率相关系数c随铁粉含量增大而减小。RG5、RG4、RG3标本相位峰值出现在高于1024Hz频率处，根据τ与特征点频率成反比关系，从相位特征点频率也能大致估算此3块标本τ值有较小数量级。3块标本相位峰值随铁粉含量的增大而增大。

含黄铜粉人工标本RG33和RG22，其m值均较大，τ值随黄铜含量的增大（3%～8%）而增大（增大一个数量级），c值随黄铜粉含量的增大而减小。两标本的相位峰值在测试频率范围内未出现。

含石墨粉标本的m值随石墨含量的增大而增大；τ值比含铁粉和黄铜粉标本的大出5～7个数量级，与含量相当的含铁块标本的τ值有相同数量级；c值随石墨含量增大而变小，且c值比含铁粉和黄铜粉标本的相对较小，本实验所有含石墨标本的c值均小于0.4。含石墨标本的相位峰值出现在0.01～0.1Hz频率之间，且相位峰值随石墨含量的增大有增大的趋势。

从浸染状标本的Cole-Cole模型参数反演结果看出：所有标本的m值均较大，且随导电矿物含量的增大而增大；含铁粉和黄铜粉的标本的τ值均较小，且随导电矿物含量的增大而增大；含石墨的标本的τ值均较大，但随导电矿物含量的增大，τ值的变化规律不明显，其c值较小；浸染状标本的c值均随导电矿物含量的增大而减小。

（2）颗粒状标本。颗粒状标本模型简图如图13（b）所示，黄铜颗粒呈圆柱形。颗粒状标本Cole-Cole模型参数随导电颗粒含量变化而变化的规律明显，随含量增大，ρ₀减小，m增大，τ增大，c增大。

相位峰值在测试频率范围内出现，且相位峰值随黄铜含量增大而增大。在相位主峰值左侧还出现另一个微弱的小峰值，如图16所示。RG8、RG9、RG24、RG25标本也有类似的相位曲线，在浸染状标本里没有出现两个相位峰值的现象，这种现象可能和导电矿物的颗粒大小有关。

图16　RG26和RG27标本的相位频谱曲线

（3）片状标本。片状标本模型如图13（c）所示。在频率范围内，复电阻率振幅频谱曲线的高频和低频部分分别逐渐趋于某一值，与标准的Cole-Cole模型曲线相似，这说明片状标本的频谱特征在实验频率范围内能得到较好反映，如图17所示。

图17　片状RG28标本复电阻率振幅及相位曲线

片状标本的相位峰值均出现在测试频率范围内，且峰值相位随含量的增大而增大。RG28和RG29为含薄片标本，连通较好，其相位曲线只有一个峰值；RG30标本［图13（e）］黄铜含量与RG28相同，但其内黄铜片较细、成网状分布。从图18看出，RG30相位曲线与颗粒状标本类似，有两个相位峰值。

图18　片状标本与网状标本相位频谱曲线

2.导电矿物连通程度与标本Cole-Cole模型参数的关系

人工标本连通程度对Cole-Cole模型参数的影响较明显，如标本RG22、RG26和RG31为含量相当的含黄铜类标本，其导电矿物连通程度从差到好，τ值能增大几个数量级。含铁类标本RG5和RG13也反映了与含黄铜类标本相同的规律，如图19所示。

3.导电矿物颗粒大小与标本Cole-Cole模型参数的关系

矿物颗粒大小是影响岩石、矿石Cole-Cole参数的重要因素，已有不少文献和文章^［1～4］讨论过τ值大小与矿物颗粒大小的关系，得出的结论是τ值的大小主要受矿物颗粒大小和连通的程度控制的。本次实验的结果也表明τ值受颗粒大小和连通程度影响较为明显。

图19　人工标本τ值与导电矿物连通程度的关系

本实验对含量相同、颗粒大小不同的含黄铜颗粒标本和含不锈钢金属球标本进行了复电阻率测试。随颗粒尺寸的增大，含黄铜颗粒标本零频电阻率ρ₀变化规律不明显；τ值增大；c值随颗粒的增大有下降的趋势；极化率m和相位峰值随颗粒增大的变化规律不明显。

含不锈钢金属颗粒的标本（RG8与RG9），各参数随颗粒半径增大而变化规律与含黄铜颗粒标本相同。

3.4　小结

（1）各类标本的Cole-Cole模型参数随其导电矿物的成分、结构、含量和颗粒大小等因素变化表现出一定的变化规律：含矿与不含矿标本的Cole-Cole模型参数有显著的区别，主要体现在ρ0、m和τ3个参数上；石墨标本的τ值和c值与浸染状标本的τ值和c值有较明显区别，本实验根据τ值和c值可以区分浸染状标本与含石墨标本；导电矿物为浸染状的标本与导电矿物为块状的标本的极化率m均较大，但相同条件下后者的τ值能比前者的大3～6个数量级，表明τ值能较好地反映标本导电矿物的连通程度。在4个Cole-Cole模型参数中，时间常数τ提供了导电矿物成分、结构和颗粒大小等信息。

（2）导电矿物为中等颗粒大小或网状分布的人工标本在频率范围内出现两个明显的相位峰值，而其他类型的标本未出现此特征，表明在一定的频率范围内有可能根据相位频谱曲线的不对称性区分大颗粒尺寸的标本。

（3）岩石、矿石的复电阻率振幅频谱和相位频谱含有丰富的信息（包括岩石、矿石的金属含量、成分和结构等），深入研究岩石、矿石的复电阻率振幅频谱和相位频谱对电法勘探的发展有较大的理论意义和实用价值。

参考文献

［1］Vanhala H，Peltoniemi M.Spectral IP studies of finish ore prospects［J］.Geophysics，1992，57（12）：1 545-1 555.

［2］Seigel H O，Vanhala H，Sheard S N.Some case histories of source discrimination using time-domain spectral IP［J］.Geophysics，1997，62（5）：1 394-1 408.

［3］Pelton W H，Ward S H，Hallof P G，et al.Mineral discrimination and removal of inductive coupling with multifrequency IP［J］.Geophysics，1978，43（3）：588-609.

［4］Johnson I M.Spectral induced polarization parameters as determined through time-domain measurements［J］.Geophysics，1984，49（11）：1 993-2 003.

［5］张塞珍，顾燕庭.地球物理勘探专辑（第3辑）——电法勘探（激发极化法）［M］.北京：地质出版社.1960.

［6］吴孝国.柯尔-柯尔模型激电时间谱的正反演方法及应用［J］.地质与勘探，1989，25（9）：37-41.

［7］罗润林，李亚南.基于时域激电数据的Cole-Cole模型频谱参数反演［C］.第九届中国国际地球电磁学讨论会论文集.2009，285-288.

［8］Kemna A，Binley A，Slater L.Crosshole IP imaging for engineering and environmental applications［J］.Geophysics，2004，69（1）：97-107.

［9］罗延钟.频谱激电法及其在近几年的发展［M］.勘查地球物理勘查地球化学文集（20），北京：地质出版社，1996.

［10］罗延钟，张桂青.频率域激电法原理［M］.北京：地质出版社，1988.

［11］武炜，杨进.时域激电谱有限元数值模拟［J］.物探与化探，2007，31（1）：55-58.

［12］刘崧.谱激电法［M］.武汉：中国地质大学出版社，1998.

［13］Coggon J H.Electromagnetic and electrical modeling by the finite element method［J］.Geophysics.1971，36（1）：132-155.

［14］罗延钟，张桂青.电子计算机在电法勘探中的应用［M］.武汉：武汉地质学院出版社，1987.

［15］孟永良，罗延钟，昌彦君.时间谱电阻率法的二维正演算法［J］.地球科学——中国地质大学学报，2000，25（6）：656-661.

［16］阮百尧，村上裕，徐世浙.激发极化数据的最小二乘二维反演方法［J］，地球科学——中国地质大学学报，1999，24（6）：619-624.

［17］刘树才，何昭友，刘志新.适合地形起伏的二维有限单元数值模拟技术［J］.物探化探计算技术，2005，27（2）：131-134.

［18］毛先进，鲍光淑.2.5维电阻率成像的新方法［J］.物探与化探，1999，23（2）：150-152.

［19］徐世浙.点电源二维电场问题中傅氏变换的波数k的选择［J］.物探化探计算技术，1988，10（3）：235-239.

［20］王长春.起伏地形2.5维频谱激电正反演研究［D］.武汉：中国地质大学硕士论文，2011.

［21］石砚斌.起伏地形2.5维时间谱激电正反演研究［D］.武汉：中国地质大学硕士论文，2011.

［22］Shi Yanbin，Zhang Shengye，Meng Yongliang.2.5-D forward numerical simulation of IP method in time domain under the condition of undulate terrain［J］.International conference on Environmental and Engi-neering Geophysics，14-19June 2010.Near-surface Geophysics and Geohazards，2010，1：544-54.

［23］曾歌明.时间域谱激电法的实验研究［D］.武汉：中国地质大学硕士论文，2008.

［24］黄理善.多相岩矿石模型的复电阻率物理模拟研究［D］.武汉：中国地质大学硕士论文，2009.

［25］Huang Lishan，Zhang Shengye，Jing Rongzhong.Study of the complex resistivity of rocks and ores model using physical simulation［G］.International conference on Environmental and Engineering Geophysics，14-19 June 2010，Near-surface Geophysics and Geohazards，2010，1：513-519.

［26］Vanhala H，Soininen H，Kukkonen I.Detecting organic chemical contaminants by spectral-induced polarization method in glacial till environment［J］.Geophysics，1992，57（8）：1 014-1 017.

［27］Tripp A C，Hohmann G W，Swift C M Jr.Two dimensional resistivity inversion［J］.Geophysics，1984，49（10）：1 708-1 717.

［28］阮百尧，村上裕，徐世浙.电阻率/激发极化率数据的二维反演程序［J］.物探化探计算技术，1999，21（2）：116-125.

［29］阮百尧.视电阻率对模型电阻率的偏导数矩阵计算方法［J］.地质与勘探.2001，37（6）：39-41.

［30］万乐.我国南方岩溶石山地区地下水勘查的新方法新技术［D］.武汉：中国地质大学博士论文，2000.

［31］黄俊革.三维电阻率/极化率有限元正演模拟与反演成像［D］.长沙：中南大学博士论文，2003.

［32］Seigel H O.Mathematical formulation and type curves for induced Polarization［J］.Geophysics.1959，24（3）.547-565.