电法报告-直流电测深正演曲线

多层水平地层地电断面电测深曲线的正
演的读书报告
姓名：***
班级：061084-27
学号：***********
指导老师：***
日期：二〇一一年五月
前言 (2)
目的 (2)
任务要求 (2)
工作过程 (2)
成果 (2)
原理 (3)
§1-多层水平地层上的对称四极电测深视电阻率表示式 (3)
1.多层水平地层地面点电流源的电场 (3)
2.多层水平地层上电测深的ρs表示式和电阻率转换函数 (5)
3.电阻率转换函数的递推公式 (6)
§2-水平地层上视电阻率的滤波算法 (6)
§3-多层水平地层的电测深曲线类型 (9)
A 二层情况 (9)
B三层情况 (9)
C四层及多层情况 (9)
编程 (10)
感想 (18)
关于多层水平地层地电断面电测深曲线的正演的读书报告
前言
目的：熟悉并掌握多层水平地层地电断面直流电测深曲线的正演
任务要求：
编制适用于n层地电断面的正演电测深程序（编程环境不限制，可用C 语言，C++，VC，VB，matlab，推荐用matlab）。

每个同学计算两个标准地电模型的正演计算
第一个模型：二层G型地电模型
第一层地层电阻率10欧姆米，第一层厚度10米；
第二层地层电阻率100欧姆米
第二个模型：三层H型地电模型
第一层地层电阻率：班号（4）×100欧姆米，第一层厚度15米；
第二层地层电阻率：序号（27）×1 欧姆米，第二层厚度20米；
第三层地层电阻率：1000欧姆米
AB/2为13个：2, 3, 4.5, 6, 9, 12, 15, 20, 30, 45, 60, 90, 120 （米）。

工作过程：
先进行原理分析，再用matlab进行编程，最后小结。

成果：
用matlab实现了n层地电断面的直流电测深正演。

原理
§-1多层水平地层上的对称四极电测深视电阻率表示式
1.多层水平地层地面点电流源的电场
如图所示，水平地面下有
n 层水平地层，
各层电阻率分别为ρ1、ρ2 … ρn ； 各层厚度分别为h 1、h 2…h n-1； 各层底面到地表的距离分别为H 1、H 2…H n-1，H n →∞。

求解思路：
U → E →ρs →Ti(λr) ① 用分离变量法求方程 ② 边界条件 ③ 通解
④ 各层的电位表达式 ⑤ 利用边界条件求待定系数 A 2、A 3，…,A n B 1、B 2，…,B n-1
设地面点电流源A 的强度为I 。

为求各层 中的电位表达式，将柱坐标系的原点设在A 点，Z 轴垂直向下。

在所设条件下，电位与角无关，满足如下形式的拉普拉斯方程：
022122=∂∂+∂∂⋅+∂∂z U
r U r r U
边界条件： （1）0|221=∞
→+r z U
（2）
0|01
=∂∂=z z
U （3） R
C R I U r z 1
10
12|
2
2=≈
→+πρ （22r z R +=）
（4）I
I
H z i H z i U U =+==||111
（5）
I I H z i i H z i z
U z U =++=∂∂=∂∂|1|11
111i ρρ
由分离变量)()(),(z Z r R z r U = 得到零阶贝赛尔方程，其通解为：
[]
⎰∞
-+=001)()()(dm
mr J e m B e m A U mz mz
式中：A (m), B(m) 为待定的积分变量m 的函数；J 0(mr) 为零阶贝赛尔函数。

利用边界条件，可以得到第一层电位公式：
⎰∞--⎥⎦
⎤⎢
⎣⎡++=0
0111)())((22dm mr J e e m B e I U mz mz mz π
ρ 第二层以下至n-1层，第 i 层的电位为： []
⎰∞-+=0
0)()()(dm mr J e m B e m A U mz i mz i i
第n 层内的电位表达式，由 0
z U
→∞
→ 得：
()100
()()n mz n U A m e J mr dm
∞
-=⎰
电测深只在地面工作，即z ＝0，故只需求出B1，
11100
2()()2I U B mr J mr dm
ρπ∞⎡⎤
=+⎢⎥⎣⎦⎰
式中：J ０(mr)为零阶第一类贝赛尔函数；B １(m)为积分变量m 的函数。

对于层数确定的水平地层，根据地层界面上电位和电流密度法向分量连续的边界条件，可具体求出B 1(m)的表示式。

例如，最简单的二层水平地层，利用ρ1 和ρ2 岩层分界面的相应边界条件可具体求出 1
1
2122121)2(1
12)(mh mh e K e K I m B ---⋅
=πρ ①
n=3时，
1121122
22((3)12231
22()212231223)
()1mh m h h mh m h h mh K e K e B m K e K e K K e --+--+-+=--+ 其中2
32323ρρρρ+-=K
2.多层水平地层上电测深的ρs 表示式和电阻率转换函数
令 )(2)(11m B I m B π
ρ= ②
则地面上电位公式为： []dm mr J
m B I U
)()(2120
1⎰∞+=π
ρ
若采用MN →0的装置测量，相应的ρs 表达式为： ⎪
⎭⎫
⎝⎛∂∂-==r U I r I E r S 22
22ππρ
[]mdm
mr J m B r S )()(2110
21⎰∞
+=ρρ
令 )](21[)(11m B m T +=ρ ③ 则多层水平地层上的电测深ρs 公式简写成：
mdm
mr J m T r S )()(1012⎰∞
=ρ
式中，T １(m)定义为电阻率转换函数又称核函数。

可见，电阻率转换函数与各层的层参数(厚度和电阻率)及积分量m 有关。

3.电阻率转换函数的递推公式
对于二层水平地层情况，若将①式先后代入②式和③式，便得到二层水
平地层的电阻率转换函数:1
1
2122121
)2(1
11)(mh mh e
K e K m T
---+=ρ 归纳每一层的电阻率转换函数，就可导出电阻率转换函数的递推公式：
=-++++-=-+--+-)
1)(()1()1)(()1()(212212i
i
i
i
mh i mh i mh i mh i i i e m T e e m T e m T ρρρ1i
1()()
()1()()/i i i i i i th mh T m T m T m th mh ρρ+++=+ ④ n
n m T ρ=)(
式中 1
1
)(22+-=i i mh mh i e e mh th 为双曲函数
i 1()()
()1()()/i i i i i i T m th mh T m T m th mh ρρ+-=
-
由此可知：当给定n ，ρ１，ρ２．．．ρn-1,ρn 和h 1，h 2，．．．h n-1。

可递推出上层或下层的电阻率转换函数。

电阻率转换函数递推公式④的导出，免去应用边界条件解方程组求系数B １(m)的计算。

§-2 水平地层上视电阻率的滤波算法
在电测深视电阻率的表达式中的被积函数可以分为两部分：一是包含地
下各层所有信息（厚度及电阻率）的电阻率转换函数T 1（m ）；此外是与地层参数无关的贝塞尔函数。

其虽没有解析计算结果，但可用线性滤波方法进行计算。

根据电阻率转换函数的变化规律，对m 的抽样取对数等间隔比较合适，因此，首先令e x =mr ，则电阻率的表达式变为
⎰∞
∞-=dx e J e r
e T r x x
x s )()()(121ρ
根据采样定理，一个函数可以用它在等间距离散抽样点上的抽样值表达： ∑
∞
-∞
=∆∆-∆∆-∆=
k x
x k x x x k x x k f x f /)(]/)(sin[)
()(ππ 将电阻率装换函数用它在x 数轴上的离散抽样值表达为：
∑∞
-∞
=∆∆∆-∆∆-=k x k x x x k x x x k x r e T r e T /)(]
/)(sin[)()(11ππ
记 ⎰∞
∞
-∆∆-∆∆-=dx e J e x
x k x x x k x C x x
k )(/)(]/)(sin[12ππ
则 ∑∞
-∞
=∆=k k x
k s C r e T r )()(1ρ
将C k =预先计算出来，实际上取有限项求和就可以达到足够的精度了。

从频谱分析的观点看，当电阻率装换函数用它在x 数轴上的离散抽样值表达式时，相当于滤取了频率高于1/2Δx 的谐波成分，因此这种计算视电阻率的方法称为滤波计算方法，C k 称为滤波系数。

为了提高线性滤波计算的精度，减少滤波系数的数目，需要对x 的抽样点进行位移，实际使用的线性滤波计算视电阻率的公式为
∑=+∆=2
1
)()(1k k k k s
x k s C r e T r ρ
式中，ρs （r ）为供电极距为r 时的视电阻率；T 1(e k Δx+s /r)为m k =e k Δx+s /r 时的电阻率转换函数；C k 为第k 个滤波系数；Δx 为抽样间隔；s 为位移系数。

实际线性滤波计算常用的抽样间隔有两种。

一种为六点是抽样间隔，即Δx=（ln10）/6=101/6，直直流电阻率测深曲线一般都采用这种抽样间隔进行计算。

另一种为十点式抽样间隔，即 Δx=（ln10）/10=101/10，电磁测深曲线一般都采用这种抽样间隔进行计算。

下表为常用的一套六点式抽样间隔的滤波系数，共有20个系数：k=1~20，其位移系数s=-2.1719，e-2.1719=0.11396。

用上述公式就可计算某个供电极距r 的视电阻率，只要计算20个对应于这个供电极距r的不同m值的电阻率转换函数T1（m k），将其与下表中对应的20个滤波系数相乘再求和就可以了。

因为要计算一条视电阻率测深曲线就需要计算多个供电极距r的视电阻率，为了减少计算工作量，取r i=e iΔx=10k/6，对应的m k=e kΔx+s/r i=e(k-i)Δx+s，这样，计算不同供电极距的视电阻率所需要的电阻率转换函数大多数可以共用。

下面是用数值滤波法计算视电阻率测深曲线的计算机流程：
（1）输入层参数，包括层数n、各层的层厚度和电阻率；
（2）输入要计算的供电极距范围，并由此得到r=e iΔx中i的变化范围：i min ~ i max；
（3）计算k-i的变化范围：（k-i）min ~（k-i）max；
（4）用电阻率转换函数递推公式，循环计算m j=e jΔx+s，j=（k-i）min，…，（k-i）max时的T1（m j）；
（5）用滤波公式，循环计算r i=e iΔ，i=i min，…，i max是的ρs（r i）。

§-3多层水平地层的电测深曲线类型
A 二层情况
二层地点断面上层电阻率为ρ１，下层电阻率
为ρ２。

有如下两种电测深曲线类型：（如右图）
（1）D型：ρ１>ρ２，基底为低阻；
（2）G型：ρ１<ρ２，基底为高阻。

B 三层情况
三层断面共有以下4种类型：（如右图）
（1）A型：ρ1<ρ2<ρ3，电阻率递增；
（2）K型：ρ1<ρ2>ρ3，中间高阻层；
（3）H型：ρ1>ρ2<ρ3，中间低阻层；
（4）Q型：ρ1>ρ2>ρ3，电阻率递减。

C 四层及多层情况
四层及多层断面电阻率测深曲线类型用三层断面类型的组合表示。

四层断面及其电阻率测深曲线类型共有8种，例如ρ1>ρ2<ρ3<ρ4的四层断面及其电阻率测深曲线类型为HA型。

五层断面及其电阻率测深曲线类型共有16中，例如ρ1<ρ2>ρ3<ρ4>ρ5的五层断面及其电阻率测深曲线类型为KHK 型。

每多一层，曲线类型增加一倍。

N层地电断面的电阻率测深曲线类型数为2n-1。

编程
要求计算两个标准地电模型的正演：
第一个模型：二层G型地电模型
ρ1=10 Ω·m h1=10 m
ρ2=100 Ω·m
第二个模型：三层H型地电模型
ρ1=班号（4）×100=400 Ω·m h1=15 m
ρ2=序号（27）×1 =27 Ω·m h2=20 m
ρ3=1000 Ω·m
AB/2为13个：2,3,4.5,6,9,12,15,20,30,45,60,90,120
程序设计流程图
程序如下：
function [ps]=dcszhengyan(ab2,n,p,h)
c=[0.003042 -0.001198 0.01284 0.0235 0.08688 0.2374 0.6194 1.1817 0.4248 -3.4507 2.7044 -1.1324 0.393 -0.1436 0.05812 -0.0252 0.01125 -0.004978 0.002072 -0.000318];%滤波系数C，共20个
rn=length(ab2);%极距的个数
d=log(10)/6;%取样间距
s=-2.1719;%位移系数
ps=zeros(rn,1);%视电阻率
for i=1:rn
for k=1:20
m=exp(k*d+s)/ab2(i);
T2=p(n);%电阻率转换函数
for j=n:-1:2
T1=p(j-1)*(p(j-1)*tanh(m*h(j-1))+T2)/(p(j-1)+T2*tanh(m*h(j-1)));%电阻率转换函数
T2=T1;
end
ps(i)=ps(i)+T1*c(k);
end
end
loglog(ab2,ps,'-o')
axis([1 1000 10 1000])
grid on
% xx=2:120;
% yy=spline(ab2,ps,xx);%三次样条插值，使曲线变光滑
% loglog(ab2,ps,'-o',xx,yy);
% axis([1 1000 10 1000])
% grid on
xlabel('AB/2 (m)');
ylabel('ρs (Ωm)');
legend('正演结果');
title('电测深曲线');
return
在Command Window中输入：
clear
close all
clc
ab2=[2,3,4.5,6,9,12,15,20,30,45,60,90,120];%极距AB/2
n=2;%两层地电断面模型
p=[10 100];%每一层的电阻率
h=[10 0];%每一层的厚度
% n=3;%三层地电断面模型
% p=[400 27 1000];%每一层的电阻率
% h=[15 20 0];%每一层的厚度
[ps]=dcs_zhengyan(ab2,n,p,h)
fid=fopen('psdate.txt','w');% open the file with write permission fprintf(fid,'AB/2 (m) ρs (Ωm)\r\n');
for i=1:length(ab2)
fprintf(fid,'%4g %12.4f\r\n',ab2(i),ps(i));
end
fclose(fid)
二层G 型地电模型正演电测深曲线图如下： 输出结果为：
10
010
1
10
2
10
3
101
10
2
10
3
AB/2 (m)
ρs (Ωm )
电测深曲线
正演结果
三层H 型地电模型正演电测深曲线图如下： 输出结果为：
10
010
1
10
2
10
3
101
10
2
10
3
AB/2 (m)
ρs (Ωm )
电测深曲线
正演结果
下面是经过三次样条插值光滑处理后的图形：
二层G 型地电模型正演电测深曲线图
10
10
10
10
3
1010
10
AB/2 (m)
ρs (Ωm )
电测深曲线
三层H 型地电模型正演电测深曲线图
10
10
10
10
3
1010
10
AB/2 (m)
ρs (Ωm )
电测深曲线
下面是增加极距AB/2个数后的处理结果：
ab2=[1 2 3 4.5 6 9 12 15 20 30 45 60 90 120 150 200 300 450 600 900 1200 1500 2000 3000 4500 6000 9000 12000 15000 20000 30000 45000 60000 90000 120000 150000];%极距AB/2
二层G 型地电模型正演电测深曲线图
10
10
10
10
10
10
5
1010
1010
AB/2 (m)
ρs (Ωm )
电测深曲线
三层H 型地电模型正演电测深曲线图
10
10
10
10
10
10
5
1010
1010
AB/2 (m)
ρs (Ωm )
电测深曲线
三层电测深曲线具有等值性：
对于H 型曲线：当ρ1、h 1和ρ3相同时，若在一定范围内按比例改变h 2和ρ2，保持S 2=h 2/ρ2值不变，则地中电场分布变化甚小，因而测点附近的电流密度j MN 也无明显改变。

这将导致不同的地电断面对应形状几乎相同的ρs 电测深曲线。

这就是H 型曲线的S 等值性。

10
10
10
10
10
10
1010
1010
AB/2 (m)
ρs (Ωm )
电测深曲线
ρ2=27×1.5=40.5 Ω·m h 2=20 ×1.5=30m
10
10
10
10
10
10
1010
1010
AB/2 (m)
ρs (Ωm )
电测深曲线
ρ2=27×2=54 Ω·m h 2=20 ×2=40m
10
10
10
10
10
10
5
1010
1010
AB/2 (m)
ρs (Ωm )
电测深曲线
ρ2=27×3=81 Ω·m h 2=20 ×3=60m
10
10
10
10
10
10
5
1010
1010
AB/2 (m)
ρs (Ωm )
电测深曲线
ρ2=27×4=108 Ω·m h 2=20 ×4=80m
感想
本次作为《电法勘探资料处理》这门课的作业，要求我们随便用一种编程工具来做多层水平地层地电断面电测深曲线的正演。

我所编写的程序是基于matlab这样的一个高层次平台，所要实现的功能是多层水平地层地电断面电测深曲线的正演模拟，绘制出电测深曲线。

因为在大二时曾经有一定的编程基础，是基于VC++平台的，后来还获得了全国计算机等级考试（二级C）优秀。

再后来到大三上，我们又学过数值分析及信号分析与处理等课程，并多次涉及到matlab的相关知识，最后我还报名参加了数理学院开设的“计算方法数学实验”自主学习课程，主要就是针对matlab的各种应用。

对比以前编程的经历，有些感想和收获。

最深的印象是要编好专业软件需要先研究好专业知识，相对来说编程实现反而比较容易。

比如我们在构建电测深曲线时，必须做到生成的电测深曲线符合地质规律，利用电阻率转换函数得到的视电阻率序列也要符合地质规律。

为了使得到的电测深曲线更近完美，就必须做好原理分析，这里我先通过阅读程志平编写的《电法勘探教程》中水平层状大地对称四极电阻率测深曲线这一节，了解到理论上多层水平地层的视电阻率表达式的推导思路，即先求得电势U，再求得电场E，最后就可得到视电阻率ρs，而首先要求得电势U的话，就得结合边界条件用分离变量的办法解出U所满足的拉普拉斯方程，又由于电测深工作在地面上进行，只需求出地面上（z=0）的U1，因此只要能求得B1（m），便可得到U1，而不一定要求所有的通解中的待定系数。

由此引出电阻率转换函数T及其递推公式的推导。

最后ρs就可以简单地表示
成由电阻率转换函数T1（m）和贝塞尔函数J0（mr）积分所得。

对于视电阻率的表达式没有解析计算结果，但可用线性滤波方法进行计算。

这也是后面如何用程序实现计算过程的关键之所在，因此用matlab编程的时候直接引用线性滤波的方法进行计算就可轻而易举地实现多层水平地层地电断面电测深曲线的正演模拟。

最后在画电测深曲线的时候，鉴于视电阻率随极距的变化规律，需选取双对数坐标。

又如在构建地电模型和极距的选取上，虽然老师直接给出了要求，但是在做的过程中我明显感觉选取极距的个数之不足，以至于尾段根本没有很好的显示出来。

这些就是我通过这次的作业认识到的编写专业软件和非专业软件最大的不同，你需要花较多的时间去把专业知识弄明白。

当然这所有的前提是你要对数值计算方法很熟悉。

在程序编写方面，matlab是一种高层次的平台，有大量的函数可以直接调用，使得很少的代码就可以实现很丰富的功能，而VC是相对基础一些的编程语言，编写的代码稍微多些。

这使得我对matlab兴趣大增。

同时我也认识到matlab中我还有很多知识需要学习，如下一个作业中要用到的可视化界面设计。

现在作业是完成了，但是我收获的并不只是那一串串程序代码，还有读书的重要与经验，这些对我今后的学习和帮助都十分重要。

以上是做多层水平地层地电断面电测深曲线的正演模拟过程中的一些感想和收获。

参考文献：《电法勘探教程》——程志平
《应用地球物理学原理》——张胜业潘玉玲。