直流电测深正演报告

直流电测深正演报告
一．原理
（1）电阻率测深法原理
电阻率测深法简称电测深，是用来探明水平层状(或近水平层状)岩石在地下分布情况的一组电阻率法变种。电测深法的装置特点是保持测量电极MN的位置固定，在不断增大供电电极距的同时，逐次进行观测。
通常要求满足以下条件：
每个测点的电测深观测结果，绘制成一条视电阻率ρs随极距AB/2变化的电测深曲线。通常将电测深曲线绘在双对数坐标纸上，其横坐标表示供电极距AB/2，纵坐标表示相应的视电阻率值。电测深曲线反映了测点下方垂直方向上电性层的变化情况。
for i = (nlayer-1):(-1):1
tt1 = 1 - exp(-2*m*thick(i));
tt2 = 1 + exp(-2*m*thick(i));
T = res(i)*(res(i)*tt1+T*tt2)/(res(i)*tt2+T*tt1);
end
B(j) = T;
end
%
%计算视电阻率
%
%画图
%
loglog(rr,rhoa,'r-o',...
'LineWidth',2,...
'MarkerEdgeColor','k',...
'MarkerFaceColor','b',...
'MarkerSize',10)
grid on;
axis([1 1000 1 1000]);
xlabel('AB/2 / m');ylabel('RHO / ohm');
'LineWidth',2,...
'MarkerEdgeColor','k',...
'MarkerFaceColor','b',...
'MarkerSize',10)
grid on;
axis([1 1000 1 1000]);
xlabel('AB/2 / m');ylabel('RHO / ohm');
function [rho] = sdcs1dford(ab2,nlayer,res,thick)
% sdcs1dford
%电测深一维正演程序
%输入参数：ab2--AB/2(m)
%
%
% c--系数，来自。。。
%
c=[0.003042,-0.001198,0.01284,0.0235,0.08688,0.2374,0.6194,1.1817,0.4248,-3.4507...
for i = (nlayer-1):(-1):1
tt1 = 1 - exp(-2*m*thick(i));
tt2 = 1 + exp(-2*m*thick(i));
T = res(i)*(res(i)*tt1+T*tt2)/(res(i)*tt2+T*tt1);
end
B(j) = T;
end
%
%计算视电阻率
求解思路：
U，E，ρs，Ti（λr）
①用分离变量法求方程
②边界条件
③通解
④各层的电位表达式
⑤利用边界条件求待定系数
a2、a3，…,anb1、b2，…,bn-1
电测深只在地面工作，即z＝0，故只需求出b1，
2、视电阻率的表达式
式中：J０（ｍｒ)为零阶第一类贝赛尔函数；B１（ｍ）为积分变量m的函数
对于层数确定的水平地层，根据地层界面上电位和电流密度法向分量连续的边界条件，可具体求出B1（m）的表示式。例如，最简单的二层水平地层，利用ρ１和ρ2岩层分界面的相应边界条件可具体求出
end
%
%输出数据到文件gmod1obs.dat
%
fid = fopen('gmod1obs.dat','w');
for iab = 1:length(rr)
fprintf(fid,'%12.5g %12.6g\n',rr(iab),rhoa(iab));
end
fclose(fid)
%
%画图
%
loglog(rr,rhoa,'r-o',...
%
rho = 0;
for j=1:20
rho = rho + B(j)*c(j);
end
return
三．成果
直流电测深正演的二层G曲线：
直流电测深正演的三层H曲线
res(1) = 10;
res(2) = 100;
thick(1) = 10;
thick(2) = 0;
rr = [
2
3
4.5
6
9
12
15
20
30
45
60
90
120
];
nab = length(rr);
for iab = 1:nab
[rho] = sdcs1dford(rr(iab),nlayer,res,thick);
2.7044,-1.1324,0.393,-0.1436,0.05812,-0.02521,0.01125,-0.004978,0.002072,-0.000318];
%
%计算T函数
%
for j = 1:20
T = res(nlayer);
m= exp(j*log(10)/6-2.1719)/ab2;
thick(2) =20;
thick(3) = 0;
rr = [
2
3
4.5
6
9
12
15
20
30
45
60
90
120
];
nab = length(rr);
for iab = 1:nab
[rho] = sdcs1dford(rr(iab),nlayer,res,thick);
rhoa(iab) = rho;
(一)水平层状大地上对称四极测深
1．多层水平地层上的对称四极电测深视电阻率表示式
(1)多层水平地层地面点电流源的电场
如图6.1-39所示，水平地面下有n层水平地层，各层电阻率分别为ρ１、ρ2…ρｎ；各层厚度分别为h１、ｈ２…ｈｎ－１；各层底面到地表的距离分别为H１、Ｈ２…Ｈｎ-１，Ｈn→∞。
设地面点电流源A的强度为I。为求各层中的电位表达式，将柱坐标系的原点设在A点，Z轴垂直向下。在所设条件下，电位与角无关，满足如下形式的拉普拉斯方程：
2.7044,-1.1324,0.393,-0.1436,0.05812,-0.02521,0.01125,-0.004978,0.002072,-0.000318];
%
%计算T函数
%
for j = 1:20
T = res(nlayer);
m= exp(j*log(10)/6-2.1719)/ab2;
归纳每一层的电阻率转换函数，就可导出电阻率转换函数的递推公式：
电阻率转换函数递推公式(6.1-66)的导出，免去应用边界条件解方程组求系数B１(m)的计算，开辟了正演计算层状大地电测深曲线的新领域。
二．Matlab程序
1.二层G型曲线
clear
close all
clc
format short g
nlayer = 2;
(2)多层水平地层上电测深的ρs表示式和电阻率转换函数
若采用MN→０的装置测量，相应的ρs表达式为：
令
则多层水平地层上的电测深ρs公式简写成：
式中，T１（ｍ）定义为电阻率转换函数又称核函数。
•可见，电阻率转换函数与各层的层参数(厚度和电阻率)及积分量m有关。
电阻率转换函数的递推公式
对于二层水平地层情况，若将(6.1-58)式先后代入(6.1-59)式和(6.1-62)式，便得到二层水平地层的电阻率转换函数:
function [rho] = sdcs1dford(ab2,nlayer,res,thick)
% sdcs1dford
%电测深一维正演程序
%输入参数：ab2--AB/2(m)
%
%
% c--系数，来自。。。
%
c=[0.003042,-0.001198,0.01284,0.0235,0.08688,0.2374,0.6194,1.1817,0.4248,-3.4507...
rhoa(iab) = rho;
end
%
%输出数据到文件gmod1obs.dat
%
fid = fopen('gmod1obs.dat','w');
for iab = 1:length(rr)
fprintf(fid,'%12.5g %12.6g\n',rr(iab),rhoa(iab));
end
fclose(fid)
%
rho = 0;
for j=1:20
rho = rho + B(j)*c(j);
end
return
2、三层H型曲线
cleat short g
nlayer = 3;
res(1) = 100;
res(2) =26;
res(3) = 1000;
thick(1) = 15;