椭球体上双频激电法的正演与反演算法

强建科，何继善

(中南大学 信息物理工程学院，湖南 长沙，410083)

摘  要：基于柯尔-柯尔频散理论，采用中间梯度上椭球体的复电阻率计算相应的幅频率，讨论椭球体的充电率、中心埋深和长轴倾角对幅频率异常和曲线形态的影响；研究幅频率资料的反演算法。研究结果表明：由于中间梯度剖面数据量非常有限，选择最小二乘迭代拟合算法是合适的；偏导数矩阵采用差分的方法求得，对模型参数进行无量纲处理，提高了反演效率；该反演算法迭代速度快，能够稳定收敛，反演效果较好。

关键词：激电法；双频激电；椭球体；幅频率反演；柯尔-柯尔频散

中图分类号：P319         文献标识码：A         文章编号：1672-7207(2007)06-1199-07

Algorithm of forward and inversion of dual-frequency induced polarization method on elliptical sphere

QIANG Jian-ke, HE Ji-shan

(School of Info-physics and Geomatics Engineering, Central South University, Changsha 410083, China)

Abstract: Based on the theory of Core-Core dispersion, the apparent amplitude frequency was computed using the complex resistivity on elliptical sphere with central gradient array. And some important factors were discussed, such as the charge-ability, the depth of center and the obliquity of long axis, which influence the apparent amplitude frequency greatly or make the form of curve distort. The inversion method of apparent amplitude frequency data was introduced. Jacobian matrix and the coefficient of linear equation were computed by the finite difference method. The results show that the least-square iterative algorithm is optimum for the apparent amplitude frequency data with the central gradient array because the profile data are limited. The model parameters normalized are taken out of their units, which can be used to quicken the inversion velocity. This inversion algorithm has some advantages, such as quick inversion speed, steady convergence and so on.

Key words: induced polarization; dual-frequency induced polarization; elliptical sphere; apparent amplitude frequency inversion; Core-Core dispersion

双频激电法^[1-3]是地面普查和详查勘探中寻找金属矿产的有效物探方法，但目前该方法在数值模拟和野外资料处理方面还存在很多不足，如对许多不规则形体还无法得到它的解析解或数值解，资料处理仅局限于一些定性解释方法和一些经验方法上，使得解释结果与实际情况差别较大。因此，急需一套有效的正演模拟和反演算法指导理论研究和野外生产实践。Agunloye^[4]研究了岩石矿物的非线性电极化的模拟和反演，Xiong^[5]研究了二层非各向同性模型中三维体的激电效应和电磁效应正演模拟算法，该算法采用积  分方程求解出频率域的极化率，与实际情况相符。Oldenburg等^[6-8]从seiger的等效电阻率出发，通过直流电阻率法的正演模拟算法，计算出无激电效应的电位和有激电效应的电位，合成视极化率，再通过迭代算法实现激电数据反演。罗延钟等^[9]从频率域激电法的幅频特性出发研究了频谱资料的最优化反演算法，阮百尧等^[10-15]依据seiger的激发极化理论，将电阻率转化成等效电阻率从而实现视极化率的二维正反演算法。强建科等^[16]研究了三维极化体的激电资料定量解释问题。目前，激电法的正演模拟和反演算法都是基于直流电阻率法正反演算法而发展起来的时间域激电数值模拟方法。对于频率域激电法，除对一些规则形体能够计算出其幅频异常外，大多数形体的异常无法了解，其反演问题更加复杂。Xiong等^{[5, 9, 12-13]}研究了一些频谱资料的正反演算法，但对双频激电法资料的处理还处于定性解释阶段。因此说，双频激电法的正反演问题是目前双频激电法进一步发展的瓶颈，深入研究该问题既有利于双频激电法的理论发展，也有利于进一步提高野外激电数据的解释水平。在此，本文作者基于柯尔-柯尔频散理论，选择椭球体模型作为研究对象，一方面改变椭球体的倾角以模拟不对称的幅频率异常，另一方面改变椭球体的长短轴以模拟一些近似板状体异常，以便模拟结果更接近实际  情况。

1  正演问题

1.1  椭球体上幅频率F_s的计算

设在地面水平均匀极化大地中，有1个旋转椭球状极化体，短轴长度为a，长轴长度为2a，相对地面的倾角为a，其中有2个坐标系(见图1)：xOy是关于椭球体本身的坐标系，可以旋转，x′Oy′是与地面平行的坐标系，它与地面测点位置联系在一起。这2种坐标系之间具有下列转换关系：

图1  旋转椭球体的坐标体系示意图

Fig. 1  Reference frame of elliptical sphere

当采用供电线AB与旋转轴共平面的中梯装置时，可以得到

这些参数的具体含义见文献[17-18]。

当围岩和椭球体在外电场的激发下产生极化时，在频率域中它们的电阻率相应为

把上式复电阻率代入椭球体极化公式即可得到

由上述复电阻率可以求出高频振幅，低频振幅，再利用双频的振幅即可求出

1.2  各参数对幅频率F_s的影响

从理论上讲，椭球体的幅频率F_s受很多参数的制约，这些参数包括：围岩的电阻率ρ₁，充电率m₁，时间常数τ₁，频率系数c₁；椭球极化体的电阻率ρ₂，充电率m₂，时间常数τ₂，频率系数c₂，以及椭球体的长轴长度a，短轴长度b，长轴倾角a，中心埋深H等。正演计算结果表明，椭球体的中心埋深H、充电率m、长轴倾角a是影响幅频率F_s的3个主要因素。下面计算幅频率随这几个参数变化时的异常情况。

基本参数设置为：椭球体长轴长度A=20 cm；短轴长度B=10 cm；中心埋深H=22 cm，长轴倾角α=90?；低频频率F_D=0.3 Hz，高频F_G=3.9 Hz。对于围岩，电阻率ρ₁=10 Ω?m，充电率m₁=0.04，时间常数τ₁=1 s，频率系数c₁=0.25；对于矿体，电阻率ρ₂=10 Ω?m，充电率m₂=0.80，时间常数τ₂=100 s，频率系数c₂=0.25。

1.2.1  椭球体中心埋深H变化对幅频率F_s的影响

幅频率F_s随椭球体中心埋深的增加，其值下降很快，曲线由陡立变得较平缓，且曲线极小值变得较不明显，见图2。

H/cm: 1—22; 2—24; 3—26; 4—30

图2  椭球体中心埋深不同时的F_s曲线

Fig. 2  F_s curves at different depths

1.2.2  椭球体充电率m₂变化对幅频率F_s的影响

幅频率F_s随椭球体的充电率m₂增加很快，曲线由平缓变得较陡立，且曲线极小值变得较明显，见   图3。

比较图2和图3可知：2组F_s曲线形态和幅值很相似，即固定椭球体充电率m₂、加深中心埋深H所得的幅频率F_s与固定椭球体的中心埋深H，减小充电率m₂所得的幅频率F_s基本相似，也就是说，中心埋深H与充电率m₂具有等效现象，这种现象常常会造成双频异常的多解问题。

m₂: 1—0.4; 2—0.6; 3—0.8

图3  椭球体中心埋深为22 cm、充电率m₂不同时的F_s曲线

Fig. 3  F_s curves at depth of 22 cm and different charge-ability m₂

1.2.3  椭球体长轴a的倾角a变化对幅频率F_s的影响

幅频率F_s随椭球体长轴的倾角a变化，其曲线幅值和曲线形态发生很大变化，曲线变化较缓的一侧与椭球体长轴近似平行，在接近地面的一端往往出现明显的极小值点，见图4。

α/(?): 1—0; 2—20; 3—40; 4—70; 5—90

图4  不同倾角时椭球体的F_s异常

Fig. 4  F_s curves at different obliquities

1.2.4  其他参数对幅频率F_s的影响

电阻率ρ、时间常数τ、频率系数c等对幅频率F_s都有影响，但相对来说，这些参数对F_s的贡献在一定条件下比较小，对于识别矿产属性很有意义。值得指出的是，模型实验和研究结果^[19-20]表明，面极化时，频率系数c为0.5~1.0，一般取0.8，体极化时频率系数c为0.1~0.5，一般取0.25。频率系数c还受矿体埋深的影响很大^[1]，在同样条件下，矿体越深，c值越小；而埋藏越浅，或矿体规模越大，c值越大。需注意    的是：

a. 由于计算椭球体异常时采用双坐标系统，原点设在椭球体中心，当椭球体长轴变化时，只能计算椭球体中心埋深，不能确定顶埋深。

b. 计算表明，基于柯尔-柯尔模型的双频幅频率随着时间常数τ增大而减小，但减小的幅值不太大，因此，在反演前根据一些已知资料估计τ。但是，时间常数τ很重要，它反映了激发极化过程的张弛特性，在野外现场，矿体颗粒的紧密程度和规模与τ密切相关。就矿物本身而言，石墨的τ最大，浸染状矿石的τ比致密块状的τ要小很多。

c. c₁和c₂选为0.25是因为国内大多数金属矿都是浸染状矿体，属体极化方式。

2  幅频率的反演算法

2.1  反演原理

假设实测的幅频率F_s能够用1个体极化的椭球体模型所产生的理论幅频率拟合得出，这个椭球体模型的各项参数就近似为所求解的答案。幅频率资料的最优化反演实质上是自动进行的“选择法”，其基本思路是预先给定一组模型参数的初值X，并按椭球体的幅频率计算公式计算理论幅频率F_s(f，X)，看其与实测幅频率曲线F_s(f)的拟合情况：若这2条曲线拟合得好，则将此理论曲线的模型参数就当作大地实际参数的解释结果；若拟合误差较大，则按一定方法计算模型参数的修改量DX=(Dx₁，Dx₂，…，Dx_n)，并用新的模型参数(X+DX)重新计算理论幅频率，……。这样，反复修改模型参数X，直到理论曲线与实测幅频率曲线达到最小二乘意义下的最佳拟合(拟合差最小)为止，并以此时的模型参数作为大地的实际参数。

依据上述分析，设幅频率最小二乘的目标函数为

2.2  模型参数处理

椭球体上的幅频率与许多参数有关，而且参数的单位各异，取值范围很大，直接影响迭代反演的速度和效果，有时甚至使方程发散，得不到有效解。阮百尧^[10]对模型参数作无量纲化处理后取得了良好的效果，这里也引用该方法对模型参数进行处理。

把式(9)代入式(7)，两边同除以F_ci(x₁, x₂, …, x_j, …, x_n)，得：

式(11)为线性方程组，解此方程组可得出新的模型参数。

2.3  反演算法流程图

整个反演步骤如下(见图5)：

a. 以文件形式读取原始数据和初始模型。

b. 正演计算椭球体异常；计算实测数据与理论数据的误差、参数无量纲处理。

c. 判断是否满足误差要求？若满足则结束，否则计算偏导数矩阵。

d. 求解偏导数的广义逆矩阵。

e. 解出参数修改量DX。

f. 按公式修改参数：。

g. 返回到步骤b.，再次计算新模型的幅频率，如此反复迭代，直到满足误差要求或给定次数为止。

2.4  算  例

为了验证中梯装置上椭球体幅频率的反演算法，用正演方法给出1条中梯剖面曲线，其模型参数如表1所示。据这条幅频率曲线可求出其倾角。反演前给理论数据增加5%的随机误差。为了提高迭代速度和精度，一些对F_s贡献小的参数就直接被假定或估计。

图5  幅频率反演流程图

Fig. 5  Process of F_s inversion

表1  模型参数反演前后结果对照表

Table 1  Comparison of models parameters before and after inversion

精度采用均方误差和修改量同时控制，当均方误差小于5%或最大参数修改量小于0.005时，迭代结束。本算例经过12次迭代，反演结果见表1，可见，反演相对误差较低，满足精度要求。

3  结  论

a. 借助椭球体模型的正演公式，可以近似模拟地下不同几何尺寸极化矿体的双频激电异常；改变围岩和极化体的电性参数，认识不同极化条件下的幅频率异常特征。

b. 极化体中心埋深H和充电率m₂具有等效现象，即固定椭球体充电率m₂，加深中心埋深H所得的幅频率F_s与固定椭球体的中心埋深H，减小充电率m₂所得的幅频率F_s的形态和幅值基本相似，这常常造成双频异常解释的多解问题。

c. 采用最小二乘优化反演算法，可以根据中间梯度的幅频率F_s数据迭代得出矿体的几何参数和部分极化参数，如椭球体的长短轴长度、中心埋深、椭球体的倾角、极化率等，为进一步钻探提供理论依据。

参考文献：

[1] 何继善. 双频激电法[M]. 北京: 高等教育出版社, 2006.

HE Ji-shan. Dual-frequency induced polarization[M]. Beijing: Higher Education Press, 2006.

[2] HE Ji-shan, BAO Guang-shu. On induced polarization surveying system[J]. JCAIMM, 1986, 6(4): 1-9.

[3] 陈蜀雁, 刘永明. 双频激电法在阿勒泰山区快速找矿评价中的应用[J]. 中南大学学报: 自然科学版, 2006, 37(3): 588-592.

CHEN Shu-yan, LIU Yong-ming. Application of dual-frequency induced polarization method[J]. Journal of Central South University: Science and Technology, 2006, 37(3): 588-592.

[4] Agunloye O. Induced polarization: Simulation and inversion of nonlinear mineral electrodes[J]. Pure and Applied Geophysics, 1983, 121(1): 63-90.

[5] XIONG Zong-hou. Induced-polarization and electromagnetic modeling of a three- dimensional body buried in a two-layer anisotropic earth[J]. Geophysics, 1986, 51(12): 2235-2246.

[6] Oldenburg D W, Li Y. Estimating depth in d_c resistivity and IP surveys[J]. Geophysics, 1999, 64(2): 403-416.

[7] Li Y, Oldenburg D W. 3-D inversion of induced polarization data[J]. Geophysics, 2000, 65(5): 1931-1945.

[8] Weller A. Three-dimensional inversion of induced polarization data from simulated waste[J]. Journal of Applied Geophysics, 2000, 44(2/3): 67-83.

[9] 罗延钟, 方 胜. 视复电阻率频谱的一种近似反演方法[J]. 地球科学: 中国地质大学学报, 1986, 11(1): 93-102.

LUO Yan-zhong, FANG Sheng. An approximate inversion of the apparent complex resistivity spectrum[J]. Earth Science: Journal of China University Geosciences, 1986, 11(1): 93-102.

[10] 阮百尧. 电阻率/激发极化法测深数据的一维最优化反演方法[J]. 桂林工学院学报, 1999, 19(4): 321-325.

RUAN Bai-yao. 1-D optimization inversion method for resistivity and IP sounding data[J]. Journal of Guilin Institute of Technology, 1999, 19(4): 321-325.

[11] 阮百尧, 村上裕, 徐世浙. 激发极化数据的最小二乘二维反演方法[J]. 地球科学: 中国地质大学学报, 1999, 24(6): 619-624.

RUAN Bai-yao, CUN Shang-yu, XU Shi-zhe. Least square 2-D inversion method for induced polarization data[J]. Earth Science: Journal of China University Geosciences, 1999, 24(6): 619-624.

[12] 阮百尧, 罗润林. 一种新的复电阻率频谱参数的递推反演方法[J]. 物探化探计算技术, 2003, 25(4): 298-301.

RUAN Bai-yao, LUO Rui-lin. A new recursive inversion method of the complex-resistivity spectrum[J]. Computing Techniques for Geophysical and Geochemical Exploration, 2003, 25(4): 298-301.

[13] 蔡军涛, 阮百尧, 罗润林. 层状大地视真频参数测深曲线的反演[J]. 中南大学学报: 自然科学版, 2004, 35(4): 662-666.

CAI Jun-tao, RUAN Bai-yao, LUO Run-li. Inverse calculation on apparent intrinsic induced polarization parametric curves over multi-layered earth[J]. Journal of Central South University: Science and Technology, 2004, 35(4): 662-666.

[14] 王华军, 罗延钟, 陈玉坤, 等. 河北某大型铅锌矿电法资料的计算机解释[J]. 物探与化探, 2001, 25(2): 144-151.

WANG Hua-jun, LUO Yan-zhong, CHEN Yu-kun, et al. Computer interpretation of electric prospecting data from a large-size lead-zinc deposit in Hebei Province[J]. Geophysical & Geochemical Exploration, 2001, 25(2): 144-151.

[15] 杨 进. 激电找水资料处理解释软件系统及其应用[J]. 物探与化探, 1999, 23(5): 363-375.

YANG Jin. The data-processing and interpretation software system for IP water exploration and its application[J]. Geophysical & Geochemical Exploration, 1999, 23(5): 363-375.

[16] 强建科, 罗延钟, 熊 彬. 固定点源测深激电异常研究[J]. 地球物理学进展, 2005, 20(4): 1176-1183.

QIANG Jian-ke, LUO Yan-zhong, XIONG Bin. The study of 3-D IP anomaly by fixed point source sounding[J]. Progress in Geophysics, 2005, 20(4): 1176-1183.

[17] 傅良魁. 激发极化法[M]. 北京: 地质出版社, 1982.

FU Liang-kui. Induced polarization method[M]. Beijing: Geological Press, 1982.

[18] 罗延钟, 张桂青. 频率域激电法原理[M]. 北京: 地质出版社, 1988.

LUO Yan-zhong, ZHANG Gui-qing. Principle of frequency-domain IP[M]. Beijing: Geological Press, 1988.

[19] 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.

[20] Wong J. An electrochemical model of the IP phenomenon in disseminated sulfide ores[J]. Geophysics, 1979, 44(7): 988-989.

收稿日期：2007-03-11；修回日期：2007-05-10

基金项目：国家自然科学基金资助项目(60672042)；中南大学博士后基金资助项目(2007年)；中国地质大学(北京)地下信息探测技术与仪器教育部重点实验室开放基金资助项目(GDL0502)

作者简介：强建科(1967-)，男，陕西岐山人，博士后，副教授，从事地球物理电磁法正、反演算法研究

通信作者：强建科，男，副教授；电话：13548985729；E-mail: qiangjianke@163.com