基于简单几何体拟合的自然电位异常反演

2018-11-05李伟林

物探化探计算技术 2018年5期

涂君, 李论, 周军, 李伟林, 王成，金杨

(1. 成都理工大学地球勘探与信息技术教育部重点实验室,成都 610059;2. 电子科技大学资源与环境学院,成都 610054)

0 引言

自然电场法(self-potential，SP)是利用岩、矿石由于电化学作用在其周围产生的自然极化电场，进行找矿、填图和解决水文地质问题的一种被动源电法勘探方法[1]。该方法近年来被广泛应用于石墨矿、硫化矿等矿产勘探当中。早在上世纪20、30年代，人们就开始着手自然电法的研究[2]。Yüngül[3]奠定了利用简单形体的响应去逼近自电异常的反演解释方法基础。对于不均匀体自电异常定量反演解释方法经过几十年的发展，已形成很多分支，但一般主要分为两大类，即复杂模型的最优化反演与简单形体模型的快速反演。

所谓复杂模型的最优化反演基于严格二、三维正演模拟，此类方法适用于任意复杂形体的自电异常。由于此类方法中的正演计算引入了有限元或有限差分法等数值计算方法，导致反演计算的效率相对较低[4]。在当前技术条件下，此类反演只能获得异常体在地面的电流值，且无法得到模型电阻率的空间分布特征[5-6]。

所谓简单形体模型的快速反演是基于单独异常体的响应分布，可用简单形体的响应近似代替的原理。该方法的合理性在于勘探获得的自电异常往往是复杂的、综合的总体异常，但该总异常是可以分解为简单异常的叠加组合[7]，因而实际工作中的自电异常可以作为简单形体来反演[8]。

简单形体模型的快速反演不同于复杂模型的最优化反演，简单形体模型的快速反演理论方法推导简单、对计算机性能需要求低、计算速度快。目前一般反演技术对原始数据利用效率不高、方法不稳定以及去噪能力不强[8-9]。针对于此，笔者提出一种基于简单形体的自电异常半自动化反演策略。

1 方法原理

设异常体的中心埋藏深度为z，极化轴与地面的夹角为θ，围岩的电阻率为ρ1，异常体的电阻率ρ2，异常体表面最大电位跃变综合参数Ur，异常形体参数q，则沿X轴方向主剖面上的电位表达式为式(1)

(1)

(2)

当异常体为球体(3D)、水平无限延伸板状体(2D)和垂直有限延伸板状体(3D)时，q值分别取1.5、1和0.5[7]。图1为球体自电异常示意图。因此简单异常体的自电异常响应可以根据公式(2)进行数值模拟计算求得。

图1 球体自电异常示意图Fig.1 A sketch showing cross-sectional view, geometries and parameters of sphere

当异常体响应对X轴方向进行求二阶导数时，有：

(3)

其中：s、Vxx分别代表求导时的步长因子和二阶导数。取x=0处的Vxx(0)值，可以将式(3)改写成式(4)。

Vxx(xi,z,s,θ)=A*B

(4)

式中：

公式(4)为计算任意点处二阶导数值计算公式。

取xi=±s处的二阶导相加并除以xi=0处的二阶导，化简得到：

(5)

(6)

式中:Vxx和s的意义同上。V(x)为各点的实测数据。

至此，得到非线性方程(5)的一般表达式，方程(5)含有未知数z(异常体中心埋深)、s(二阶导数步长因子)以及q(异常体形体参数)。解方程(5)得到在不同步长因子s下的q、z值；绘制z_q曲线，出现交点的横、纵坐标即为异常体的埋深值和形体参数值q值。

利用x=0处的V(0)带入公式(2)可以得到k的表达式：

(7)

带回公式(2)重新整理得到电位表达式如下：

(8)

式中的zc表示由前面计算得到的异常体埋深。因为z和q值已知，所以公式(8)只含极化轴和地面的夹角θ未知。利用最小二乘法原理可以得到计算夹角θ的表达式为式(9)。

(9)

式中:V(xi)为实测值;V(xi,zc,q,θ)为当前模型参数下理论值。令式(9)取极小值得到θ的表达式为式(10)。

θc=cot-1(C-D)

(10)

式中：

(11)

带入zc和q值，计算公式(11)即可以得到θ的值θc。

同理将得到的zc、q、θc值设为固定值，带入公式(2)再次利用最小二乘法得到电偶极矩k的计算表达式：

(12)

从公式(5)、公式(11)和公式(12)分别计算了异常体的中心埋深z、形体参数q、极化轴与地面的夹角θ和电偶极矩k。其基本流程见图2。

图2 算法流程Fig.2 Algorithm flow

2 理论模型计算

正演计算是已知模型空间求解数据空间，而反演是已知数据空间，求取模型空间；正演是基础，反演是目的。笔者选择一个典型模型按文中方法进行反演，验证方法的正确性。

2.1 正演模拟计算

利用正演计算对理论模型自电响应特征进行分析，为观测数据的定性以及定量解释奠定理论基础。

表1 正演模型的埋深参数

图3 不同埋深模型自电异常响应Fig.3 SP anomaly aroused by different depths

模型Z/mqθ/°k/mV·mB1100.845-1000B210145-1000B3101.545-1000

正演计算中心埋深不同模型的响应并分析，参数见表1。

图3为模型埋深不同的正演响应曲线，在模型其他参数相同的情况下，异常体的埋深增加使响应极小值幅值呈非线性递减。异常极小值对应坐标位置基本保持不变，极大值所对应的坐标位置随着埋深增加而远离原点，并且在远离中心的位置，响应趋于零。

正演计算模型形体参数不同的响应特征，见表2。

图4为模型形体参数不同下的正演响应曲线，在模型其他参数相同的情况下，异常体的形体参数增加使响应极小值幅值呈非线性递减；异常极小值对应坐标位置基本保持不变，极大值对应坐标位置随着形体参数增加而靠拢异常体中心位置。整体响应可见球体(q=0.5)响应的幅值最小。

图4 不同形体参数模型自电异常响应Fig.4 SP anomaly aroused by different shape factor

模型Z/mqθ/°k/mV·mC1101.50-1000C2101.545-1000C3101.590-1000

图5 不同倾角模型自电异常响应Fig.5 SP anomaly aroused by different inclination angle

正演计算模型极化轴倾角不同的响应特征并分析(表3)。

从图5中可以得到，随着极化轴倾角的变化，自电异常也发生这变化[1]。具体表现为：

2)当极化体为倾斜极化(0°<θ<90°)时，其电位曲线介于水平极化和垂直极化电位曲线之间。倾斜极化电位曲线的形态随倾角的不同而不同：电位曲线的极小值已不在极化体正上方，而是向极化轴倾斜的反方向移动。θ越小，移动的距离越大；零值点情况亦是随θ的减小向极化轴倾斜相反方向移动。零值点与原点的距离为:x0=z*tgθ。

由此可见，z一定时，θ角越大，零值点偏离原点越远；当θ=90°时，零值点将在无穷远处。因极化轴倾斜，在倾斜一侧出现的电位正值，是极化轴倾斜较小的标志，且随着倾角θ的减小，电位正值逐渐增大。在自然界中，由于水文地质条件关系，一般极化轴近于垂直，故在金属矿体上常观测到负电位，只在地形切割很强的地区位于陡峭的金属矿体上，有时能见到显著的电位正异常。并且从公式(2)可见，电位与电偶极矩k呈正相关关系。

2.2 理论模型反演

分别设计简单参数理论模型和一个复杂参数理论模型进行验证。根据表4中模型参数，计算正演响应，并作为初始实测数据，带入反演流程，计算出的二阶导数曲线如图6所示。

表4 模型反演参数

图6 自电异常二阶导数曲线Fig.6 Second derivative of SP anomalies

得到二阶导数值后，进行非线性方程组的求解，得到固定s下的z_q曲线(图7)。图7是在不同步长因子下，反演所得关于异常体埋深与形体参数的曲线图。图7中横坐标即为异常体形体参数，纵坐标为异常体埋深参数。

图7 模型反演z_q曲线Fig.7 The inversion z_q curves

模型Z/mqθ/°k/mV·mF14.50.955.6-1950

图8 自电异常二阶导数Fig.8 Second derivative of SP anomalies

从二阶导曲线左右不对称性，可以初步判断异常体极化轴的倾向，二阶导数极小值所在坐标原点的一侧为倾向方向。从图7得到模型埋深和形体参数q分别为10 m和1.5。将得到的z、q值带入公式(11)、公式(12)即可得到、k值分别为45°和-1000 mv.m。计算结果与设置的初始模型参数一致，初步说明该方法是可行的。

对于理论模型E1进行反演后，验证了该方法对简单模型能够进行正确的反演。设计原始模型参数如表5所示，分别得到反演过程中的二阶导数曲线图和z_q曲线图。

图9 模型反演z_q曲线Fig.9 The inversion z_q curves

图8为理论模型F1的自电异常二阶导数，图9为该相应模型的反演z-q曲线。从图9可得到zc、q值分别为4.5 m和0.9；将得到的z、q值带入公式(11)、公式(12)得到θc、k值分别为55.6°和-1950 mV.m。计算结果与原始模型参数一致，进一步说明该方法是正确可行的。

3 探测实例

3.1 方法验算

通过理论模型的数值模拟计算，验证了文中所提方法的可行性。为进一步验证该方法的实用性，选取与已知工区进行对比验证。选取Yüngül在土耳其东南部某铜矿区采集一条自电剖面，并且Bhattacharya B B[10]对该剖面进行处理，Essa K[8]也处理了该剖面。笔者选取二阶导步长因子s为3、5和7，对其进行处理。先计算二阶导数，然后得到z_q曲线交汇图，从中读取埋深和形体参数值分别为=35.903 m和q=1。表明此异常体用水平无限板状体(q=1)能较好拟合。将计算出来的z、q值带入公式(11)、公式(12)，计算得到极化轴与地面的倾角为17.8201°，电偶极矩k为-12072.8 mV·m，即得到反演参数(表6)。

表6 反演结果

计算反演所得参数的正演响应，并与实际资料和原文献结果进行对比分析。由图10可见，响应结果整体吻合很好，该方法相比较于Essa K方法更能拟合原始剖面，但也存在误差，其产生原因可能是地下异常体空间展布复杂引起的。综合整条剖面，此方法整体较好的拟合了原始剖面，说明方法对实际资料也能适用。

图10 反演响应对比Fig.10 The comparison graph of inversion results

3.2 石墨矿区探测数据反演

测区位于杨子克拉通北缘，经历了结晶基底形成，褶皱基底形成，澄江湖大陆裂谷、克拉通盆地演化，内陆盆山耦合—推覆构造五大演化阶段。区内矿产主要有铁矿、钾长石矿、霞石铝矿、石墨等。富含有机质陆源碎屑及原始生物沉积形成是碳质的主要来源[11]。

石墨矿的碳源主要为有机成因生物碳，石墨矿本身具有低阻高极化特性，而石墨矿的成矿围岩多为电阻率较高、激化率极低的岩体，因而石墨和其他岩(矿)石相比具有明显的低阻高极化特性，并且在成矿后具有稳定的层位和一定的规模。石墨矿体上方自然电位值最低，勘查区岩石与矿石之间的存在电阻率和自然电位电性差异，这为物探工作的开展提供了较为理想的地球物理条件。

选取工区一条自电剖面数据进行反演解释，该剖面含测点44个，数据采集点距40 m，对实测数据进行反演前预处理，得到反演结果如表7所示，表中角度为按李金铭[1]定义的正方向，即从右到左变化。

表7 反演结果

计算反演所得参数的正演响应，并与实际测量资料(图11)、长偏移距瞬变电磁法(LOTEM)反演结果和实际地质钻进图进行对比分析(图12)。

图11中，黑色点线为经过预处理后的测量数据，原始数据稳定性较好，异常幅值明显。而图11中的红色曲线为利用本文快速反演方法得到的反演结果，反演结果曲线更为光滑，整体趋势与原始数据拟合较好。

图11 反演响应对比Fig.11 The comparison graph of inversion results

图12分别为石墨矿体上的自电异常测量数据剖面图与反演拟合数据，利用线源长偏移距瞬变电磁得到的石墨矿体反演响应图，以及利用钻井资料得到的地质剖面成果图。从图12中钻井资料得到，石墨矿体倾角为向右大角度陡立倾斜，而且石墨矿体中心点埋深也在80 m左右，形体基本可以用垂直有限板状体来拟合。从自电剖面上可以看出，自然电位低阻特征与矿体位置对应较好。地面长偏移距瞬变电磁反演结果，电阻率总体呈“两高夹一低”的特征。自电反演得到的埋深和倾角以及形体参数与瞬变电磁反演结果和钻进资料吻合较好，所以通过图12所示的勘探实例，进一步验证了该方法。

在石墨矿体的平面边界圈定方面，目前的自然电位方法发挥了重要作用，未来的石墨矿探测工作中，自然电位方法将占据重要地位。因此，有必要着力于该技术的方法机理研究，从数值模拟、实测资料的快速成像方法等方面入手[12]，进一步优化和提升自然电位数据的解释策略，使之可以更好的为石墨矿及其他矿产勘探服务。

4 结论

从理论上推导了该方法的可行性，设计理论模型对方法的正确性进行了验证，应用于实际资料反演解释中，验证了其可靠性。因为该方法是利用二阶导数解构建非线性方程，所以推导简单，易于实现。在解倾角和电偶极距时，利用剖面全部数据，因而数据使用率高，可靠性大。相比较于复杂模型的最优化反演，此方法方便、快速，适于初步快速反演。且在用该半自动反演方法时，加入了人为对剖面的地质先验信息，从而提高了数据的可靠性。