1 引言

为了有效地消除电磁耦合和评价激电异常，Pelton等^[1]提出了频谱激电法(SIP).它利用常规电阻率法的观测方式，在若干个频率(频率范围n×0.01~n×10²Hz)下测量视复电阻率谱和视相位谱，推断目标体的真复频谱参数特性，从而提高对异常客观评价的准确性.Pelton 等^[1]认为，对均匀岩、矿石，由激电效应引起的复电阻率随频率的变化可由Cole-Cole模型描述，并通过对大量岩、矿石标本的广域频谱复电阻率测量，验证了该模型的正确性，并指出可以根据岩、矿石的Cole-Cole模型参数，按结构评价极化异常体.

由于SIP 法可观测的参数多，多参数组合的解释能为评价激电异常源的性质提供更丰富的信息，且勘探成本低，周期短，抗干扰能力强，近年已在矿产勘查、油气藏探测及环境调查等方面得到了广泛应用^[2-4].尤其在隐伏多金属矿勘探方面，复电阻率法能较好地解决矿与非矿异常的区分问题，提高了寻找隐伏矿的能力.

随着SIP 法的理论研究、仪器设备及施工技术日益成熟，与之相应的反演算法也得到了快速发展.我国学者罗延钟提出了由视谱参数解析计算极化体真谱参数的方法^[5].张桂青等^[6]提出由视谱直接反演极化体真谱的思想，并建立了稀释系数为常数及随频率变化下的反演方法.刘崧等^[7]提出了联合谱激电反演真Cole-Cole参数的思想，并实现了椭球极化体的真谱参数反演.蔡军涛等^[8]在直流电阻率法数值模拟的基础上，提出了利用最优化反演和递推反演相结合求取真频参数的反演方法^[9].Loke等^[10]同样在不考虑电磁感应的情况下，以近似反演得到的非均匀半空间作为初始模型，实现了理论数据的SIP二维反演.就以上SIP反演方法而言，在面向实际应用时存在着一些难以避免的问题.利用视谱直接拟合真谱的方法，在建立视谱和真谱的近似关系式时，难以对普遍的情况使用统一的复电阻率模型进行描述，在解释复杂构造异常体及真谱参数的空间分布上存在着一定的局限性，并且将电磁效应仅用一个基于Cole-Cole模型的简单数学表达式描述时，也很可能造成解释上的偏差.而以直流电阻率法正演为基础的SIP 反演算法，一般均建立在可忽略电磁耦合这一假设之上，而实测的SIP 数据中却一定包含了此种效应，特别是在观测频率较高时，电磁效应更为显著，这意味着电磁耦合“噪声"在此类反演中将起着更大的干扰作用.

近年来，由频率域Maxwell方程出发建立的SIP反演方法得到了广泛研究.张辉等^[11]基于体积分方程的SIP正演模拟，首先提出并完成了基于阻尼最小二乘法的SIP三维反演，李建平^[12]在此基础上进行了改进，实现了起伏地形下的SIP三维反演.AhmadGhorbani^[13]编写了包含电磁效应的SIP 一维反演代码，并验证了其正确性.徐凯军^[14]将解耦的波数域电磁场偏微分方程用有限元法实现了正演模拟，并利用阻尼最小二乘法通过对电场的拟合，研究了复杂地形的2.5维SIP 局部反演.梁盛军^[15]将共轭梯度法和阻尼最小二乘法按先后顺序应用于三维SIP反演中，并对理论模型进行了局部反演计算.可以看出，虽然以上反演研究既考虑了电磁效应和激电效应并存的情况，又能直接反演出目标体的多种复电阻率参数、位置和几何特征.但是，由于反演单元存在多个Cole-Cole模型参数，在二维和三维全区反演时，模型参数量远多于观测数据量，导致反演问题的欠定性严重.因此，该方面的研究目前也仅停滞在理论研究及小规模的模型试算阶段，尚未达到实际应用的程度.

针对频谱激电反演中存在的问题，尤其是目前最符合实际应用的2.5 维反演问题，本文提出并实现了利用多个排列的视电阻率和视相位数据的联合反演方法，在考虑电磁效应的情况下，能够同时反演出二维地质断面上所有单元的复电阻率参数，解决了SIP反演欠定性严重的问题.算法中利用最小二乘原理构建了反演目标函数，并加入了Occam 反演^[16-17]的光滑模型约束，以增加反演的稳定性.对于反演中的雅克比矩阵的计算问题，本文借助了电场的偏导数形式，推导出了灵敏度元素的解析表达式，并应用互换定理^[18-20]对其进行直接求解.最后，对安徽某地区的SIP实测数据进行了2.5 维反演成像，通过与已知钻孔资料以及和CSAMT 反演结果的对比分析，验证了该反演算法的应用效果.

2 2.5维复电阻率反演原理

SIP反演问题可归结为求解以下泛函极值的问题:

(1)

其中，x为模型的复电阻率参数向量；ε 为数据拟合误差；f为实测数据向量；F(x)为正演函数.考虑到电磁效应对观测数据的影响，正演模拟是以波数-频率域的电磁耦合方程^[21]为理论基础，将实电阻率变量用(2)式的Cole-Cole模型替换后，再用等参单元的有限元法^[22]完成了离散化处理，最终可得到由式(3)表示的正演线性方程组^[23]:

(2)

式中，ρ(iω)为复电阻率；ω 是角频率；ρ₀、m、τ、c为Cole-Cole模型参数，分别表示零频电阻率、极化率、时间常数、频率相关系数.

(3)

式中，K是对称、带状的系数矩阵；F为剖分网格节点处的波数域待求场向量；B为发射源项.当求解出空间域电场值后，利用以下近似公式就能计算出测点处视电阻率及视相位值.

(4)

式中，ρ_a 和Φ_a 分别为观测点处的视电阻率值与视相位值；K是装置系数；E为观测偶极中心处的电场值；ΔL为观测偶极长度；I为发射电流.

2.1 反演方程组的建立

为了充分利用视电阻率和视相位数据包含的地质信息，本文采用了联合反演的策略.从公式(4)中不难发现，复电阻率法的反演拟合，本质上是基于电场的振幅与相位的拟合，因此，文中给出的反演算法均针对电场的振幅及相位进行推导.考虑到四种SIP参数的量级不同，为了提高反演的稳定性，对参数进行了对数归一化处理，即令x= (ln(ρ0)，ln(m)，ln(τ)，ln(c)).遵循以上规则，以Occam 原理构建了单一排列下的反演目标函数:

(5)

式中，a和Φ 分别为观测电场的振幅向量与相位向量；A(x)和Φ(x)分别为电场振幅与相位的正演函数；W_a和W_Φ 为归一化对角阵；λ 是阻尼因子；μ 为缩放系数；R为模型的二阶粗糙度矩阵.

将(5)式在k次迭代模型x_k邻域线性展开，并极小化，可得到以下实系数反演方程组:

(6)

其中，Δd_A和Δd_φ 分别为电场振幅和相位的绝对拟合差向量；J_A和J_φ 分别为电场振幅与相位的灵敏度矩阵；Δx_k=x_k+1-x_k为模型修正向量.通过求解以上方程组，就能得到本次迭代更新的模型参数向量x_k+1，再将其作为初始模型进行下一次迭代计算，直到拟合误差收敛至设定的阈值.

当利用M个排列的数据反演时，首先求出各排列对应于(6)式中左端的灵敏度相关项P₁，P₂…，P_M，及右端项S₁，S₂…，S_M，再令其线性相加生成新的灵敏度相关项Q=P₁ +P₂ + … +P_M及右端项T=S₁ +S₂ + … +S_M，此时方程组(6)变为:

(7)

通过上式便可实现利用多个排列数据的SIP 反演计算.

2.2 灵敏度矩阵推导

灵敏度矩阵的建立是反演的重要环节，它直接决定了反演成像的时间和精度.在(6)式中可以看出，反演方程中包含了J_A和J_φ 两种类型的灵敏度矩阵，在反演过程中需要分别求取，其解析表达可借助场的灵敏度形式获得.在频率域中，将空间分布的电场表示成以下复数形式:

(8)

式中，$\left.E \right|$和$\varphi $ 分别是电场的振幅和相位函数；i为虚数单位.对(8)式两端求取x_j的偏导数(x_j是第j块反演单元上任一种Cole-Cole模型参数)，经整理后有:

(9)

再根据复合函数求导法则，便可推导出以下各复电阻率参数的灵敏度解析表达式:

(10)

在(10)式中，ω 为角频率；ρ 为复电阻率变量，与σ 互为倒数；为j单元的复电阻率对其Cole-Cole 模型参数的偏导数，可解析求取^[24].这样，灵敏度元素J_Aj和J_Φj的计算，最终便归结到电场对第j块单元复电阻率的偏导数计算问题，即求解为节约计算成本，应用互换定理，并通过伪正演对该项进行计算.鉴于互换定理在2.5维反演中应用的复杂性，限于篇幅，拟另文发表.

3 安徽某斑岩铜矿上的应用试验 3.1 矿床地质

安徽某斑岩铜矿床位于庐枞火山岩盆地西北缘，矿床的形成与中国东部燕山期岩浆侵入与喷发活动有关，是东部地区一个典型的岩浆热液型矿床.矿区地层简单，主要有第四系(Q)、白垩系杨湾组(K₁y)、志留系高家边组(S₁g)和坟头组(S₂f)，早侏罗统磨山组(J₁m)和中侏罗统罗岭组(J₂l)地层以及早白垩纪龙门院组和浮山组火山岩系.

矿区构造主要表现为褶皱和断裂两种类型.高家边组和坟头组地层组成NNE 向背斜，背斜上断裂交汇处为铜矿的有利富集部位.区内岩浆活动强烈，形成了一套以石英闪长斑岩和黑云母石英闪长斑岩为主的钙碱性系列的中酸性岩体，总体呈北东向分布.侵入于志留系和侏罗系地层中的石英闪长斑岩、黑云母石英闪长斑岩等为主要容矿岩体，总体呈北东-北北东向沿背斜核部分布^[25].

矿石矿物成分简单，主要金属矿物有黄铁矿、黄铜矿、斑铜矿、辉钼矿、磁铁矿、辉铜矿；非金属矿物除原岩蚀变矿物及交代残余矿物外，还有石英和长石等.矿石构造以浸染状、细脉状和细脉浸染状为主，其中含铜斑岩型矿石普遍为浸染状矿化叠加疏密不等的细脉状矿化.

3.2 野外数据采集

数据采集时选择了一条有钻井控制的剖面，使用的测量仪器是加拿大Phoenix地球物理公司研制生产的V8多功能电法测量系统，观测方式为偶极-偶极.共进行了16组排列的数据采集，每排列12个观测道，发射偶极长度为200 m,接收偶极长度为50m,共进行25个频率(0.0313~128 Hz)的观测.反演前，对观测数据进行了筛选，剔除了部分干扰较强的排列、频点及观测道.

3.3 反演结果与分析

利用本文提出并实现的2.5维SIP 反演算法，对该地区实测数据进行了反演成像.反演拟合情况见图 1~4.其中，图 1和图 2分别为1 Hz的视电阻率和视相位拟合断面图；图 3和图 4分别为16Hz的视电阻率和视相位拟合断面图.经对比可以发现，两个频率数据的反演拟合情况大致相当，总体拟合程度较好，仅局部存在微小的差异.但视相位较视电阻率的拟合程度稍差，这是因为视相位对SIP 参数的灵敏度较视电阻率高，同时也说明了视相位数据提供了更多的SIP异常信息，在反演中起到了更为关键的作用.

四个SIP参数的反演结果见图 5，图中白线围成的封闭区域为钻井控制的矿体范围.图 6 是该剖面的CSAMT 法反演结果.由反演结果不难看出，SIP反演的零频电阻率与CSAMT 法反演的电阻率剖面吻合较好，其电阻率结构可以反映出矿区的地层、岩体及背斜构造.主体高阻区是含矿岩体的反映，而分支的高阻体应该是岩支的表现，而低阻区域在主体上应该是志留系地层的表现.对于矿体而言，主体矿体具有较高的电阻率，这与矿床的类型和含矿体的电阻率测量结果是一致的(表 1).斑岩型矿体一般具有高的电阻率，特别是石英斑岩具有较其它斑岩更大的电阻率值.SIP 反演结果显示，矿体具有较小的时间常数，这与侵染状或细脉状矿石的时间常数较小的特征较为一致.在矿体范围，频率相关系数具有由浅至深逐渐变大的趋势，这也同该参数的物理意义一致，因为该参数主要反映矿物颗粒的均一性特征.矿体的浅部由于埋藏浅，不均匀性严重，因此频率相关系数较小，而深部压力大，形成的颗粒更均质，频率相关系数会相对较大.

反演的极化率在矿体上约为20.0%，与标本测量结果基本一致，但并没有出现人们期待的高极化异常.相反，志留系地层则表现出较高的极化率特征，这同粘土质粉砂岩标本相对较低的极化率(7.0%)相矛盾.我们认为，可能是粘土质粉砂岩标本的代表性出现了问题，志留系砂泥岩可能存在较多的碳质成分，在后期岩体侵入过程中接触变质，因此应该表现为较成矿岩体具有更高的极化率特征.就该区的斑岩型铜矿而言，不难看出:该矿体具有高电阻率，相对低极化，中等频率相关系数和较小的时间常数的特点.基于如上认识，推测反演图右侧(水平坐标2600，深度约200m)的异常可能是矿致异常，同样，左侧(水平坐标1200，深度约200 m)也可能是相同类型的矿致异常.

对于剖面上其它参数组合的异常，我们也进行了分析和推断.在水平坐标1500 m,深度约200 m处存在的异常，表现为电阻率低，极化率较低，频率相关系数高和时间常数高的特点，推断为隐伏角砾岩.浅地表的相对高阻(相对于围岩地层)，高极化率，低频率相关系数和低时间常数的异常，推断为斑岩型铜矿的外围似千枚岩化蚀变带(黄铁矿化，石英和绢云母化).

4 结论

针对复电阻率法反演在实际应用中存在的问题，特别是生产中急需的2.5维反演问题，本文提出并实现了2.5维复电阻率法反演，并进行了应用试验，得到了如下结论:

(1) 本文提出并实现的复电阻率反演算法，不对电磁效应做任何假定，是建立在频率域Maxwell方程基础上的完全意义上的拟合反演.不但能同时反演出全区的零频电阻率、极化率、频率相关系数和时间常数，而且能反演出异常体的几何参数.

(2) 算法中同时利用了视电阻率和视相位数据进行联合反演拟合，在增加数据信息的基础上，有效的提高了反演的分辨率.通过拟合断面图的对比发现，视相位对模型参数的灵敏性高于视电阻率，同时也说明了视相位数据提供了更多的SIP 异常信息，在反演中起到了更为关键的作用.

(3) 利用反演程序对已知矿区的实测数据进行了SIP 反演实验，通过与CSAMT 的反演结果、钻探资料和物性资料的对比验证，表明了该反演方法能成功地圈定矿区的地质构造及矿体，使按结构区分矿与非矿成为可能，具有良好的应用效果.

(4) 从反演实例可以看出，即使SIP的多参数反演能够提供更为丰富的异常信息，但在实际生产中，地质情况往往更加复杂，仅依靠反演结果难以实现精确的地质解释.因此，需紧密结合实际地质资料进行综合分析，才能更准确地解决对矿床的定位与预测问题.