中梯装置脉状体ρ_s异常的计算.docx

中梯装置脉状体ρs异常的计算
姓名：scnjlybd 学号：xxxxxxx专业：地球物理学

【摘要】：有一位于均匀场中，电阻率为ρ2的脉状体（椭球体），其长半轴为c、短半轴为a、埋深为h，求在围岩介质为均匀各向同性，电阻率为ρ1的条件下，主剖面上椭球体的视电阻率分布。
【关键词】：中梯装置；脉状体；视电阻率异常
Calculation of Apparent Resistivity Anomaly
of Intermediate Gradient Device

【ABSTRACT】：An ellipsoid situated in the uniform field is provided of which resistivity is ρ2, semi-major axis is c, semi-minor axis is a, and burial depth is h. Now we know the resistivity of the medium in the homogeneous and isotropic surrounding rock is ρ1, and then seek for the apparent resistivity distribution on main section of the ellipsoid.
【KEY WORDS】：intermediate gradient device; ellipsoid; apparent resistivity anomaly
1 方法原理
脉状体是野外实际条件下，常见的一种地质体。对于这类形状的矿体，我们可以通过对三轴椭球体的理论计算方法或通过模型实验的方法研究ρs异常的分布特征及规律。对寻找有一定走向长度的脉状矿体而言，当AB垂直矿体走向供电时，称为纵向中梯或常用中梯；当AB沿着（或平行）矿体走向供电时便称为横向中梯。这里主要简要说明前者。
1.1 均匀电流场中椭球体的一次场解答及ρs表达式
Fig 1均匀电流场中的椭球体
设在电阻率为ρ1的均匀各向同性无限岩石中有一电阻率为ρ2，3个半轴分别为a、b、c的椭球体。我们取直角坐标系，原点位于椭球体中心，坐标轴X,Y，Z分别与a、b、c重合。均匀外电场j0的方向与X轴一致，见图1。
在椭球坐标系中解拉普拉斯方程并引用极限条件，可求解椭球内部和外部一次场的电位表达式（何裕盛等,1978）为：
U1(2)=-(1+F1L0)j0ρ1x (1)
U1(1)=-(1+F1L)j0ρ1x (2)
式中：
L0=2πabc0∞dt(a2+t)(a2+t)(b2+t)(c2t)
L=2πabc0∞dt(a2+t)(u2+t)(b2+t)(c2+t)
L中的积分下限t0为以下方程的最大实根：
x2a2+t0+y2b2+t0+z2c2+t0=1
以上为全空间情况，当有平行于XOY平面的地面（z=h0）存在时，地面一次场电位表达式，可用（2）式中异常部分加倍的办法近似地写出：
U1=-(1+2FL)j0ρ1x （3）
沿X方向的一次场电场强度为
E1x=-∂U1∂x=(1+2F1g')j0ρ1 (4)
式中：g'=∂∂x(xL)=L+x∂L∂x ,为椭球体的几何系数。此时，L中的积分下限t0是方程
x2a2+t0+y2b2+t0+z2c2+t0=1
之最大实根，其中h0为椭球中心深度。而
x∂L∂x=-4πabcx2(a2+t0)2(a2+t0)2(b2+t0)2(c2+t0)2x2(a2+t0)2+y2(b2+t0)2+h02(c2+t0)2
由（4）式可写出地面上沿X方向的视电阻率表达式为
ρs=(1+2F1g')ρ1 （5）
对于b→∞情况，（5）式中的F1，g'各量可以简化，以简化后的各量代入（5）式，便得到走向无限、柱轴平行行地面的椭球柱体ρs表达式
ρsρ11+2μ12-1aμ12+c∙aca-cc2-a22a2+t032c2+t0∙∂t0∂x∙x+1-c2+t0a2+t0 (6)
1.2 脉状体上的中梯ρs异常
可以理解，脉状体矿体上中梯ρs异常的量值大小和分布特征，有仅与其埋藏深度和相对围岩电阻率的大小有关，而且还与脉状矿体的产状有关。
当令（6）式中的a<>c,则得相当于水平薄脉的ρs表达式。对于倾斜薄脉的情况也可导出与其相应的ρs表达式。根据这些ρs表达式即可算出各种ρs剖面曲线。考虑到理论计算和模型实验结果是一致的，为简单起见，这里主要用模型实验结果，对不同产状不同电阻率的脉状矿体上之中梯ρs异常进行讨论。
1.2.1 低阻脉状体上的中梯ρs异常
对于直立铜板来说，在适当条件下，由于均匀电流场的方向与模型走向方向垂直，故在矿顶上视电阻率的变化不明显，近于和围岩（水）电阻率值（约为20Ω·m）相等。因此，用中间梯度法寻找直立良导薄矿脉是不利的。
当矿体倾斜时，ρs则出现了明显异常。其特点是在倾斜方向上ρs明显降底，在反倾斜方向上则ρs有所升高，曲线呈不对称状。根据ρs有极小值的那一边可指明矿体的倾向，并能说明矿体为低阻体。这是因为矿体倾斜时，均匀外电流场的方向与矿体斜交，因此吸引电流的能力加强，致使倾斜方向的近地表处电流密度jMN减小（小于j0）,根据ρs=jMNj0ρMN的关系式，应有ρs<ρ1的低阻异常。
1.2.2高阻脉状体上的中梯ρs异常
在高阻脉状体上，不论产状如何均有大于背景值（ρ水）的ρs高阻异常，但以直立高阻脉上的异常为最大。这是因为对高阻脉来说，当脉的产状为直立时，均匀外电流场的方向与脉的走向方向垂直，因而排斥电流的能力最强。
当高阻脉为倾斜产状时，由于排斥电流的能力减弱，所以ρs异常变小。曲线不对称且在倾斜方向上ρs下降较快，有明显低于ρ1的极小值。依此可判断高阻脉的倾斜方向，并能说明异常源为高阻体。
当高阻脉的产状为水平时，由于排斥电流的能力最差，所以ρs异常最小。若脉的厚度很薄，甚至观测不到ρs异常因此，用中间梯度法寻找高阻直立脉最有利，而寻找高阻水平脉最不利。对高阻情况，视电阻率极大值点约位于矿顶上，两侧有次级极小点（倾向一侧较明显）；而低阻体情况则不同，矿顶附近对应的视电阻率异常为零，两侧有量值相近的极小（倾向一侧）和极大（把倾向一侧）约呈倒像对应。视电阻率异常形态随ρ2的变化是由于矿内、外异常电流的分布为装置类型和矿体的几何特征所决定的。虽然表达公式较为繁琐，但通过计算可以明确其间之相应关系，这无论对理论研究和野外资料解释均有实际意义。
2 实现过程
2.1建立数学模型
电阻率为ρ2的椭球体矿体，其长半轴为c、短半轴为a、埋深为h，其围岩介质为均匀各向同性，电阻率为ρ1如图2所示。观测点位于主剖面上的视电阻率表达式：
ρs=ρ11+2ρ1-ρ2aca-cρ2a+ρ1c1-c2+ua2+u-12a2-c2(c2+u)12a2+u23∂u∂xx (3-1)
其中
u=12(h2+x2-a2-c2)+12x2+h22+a2-c22-2(a2-c2)(x2-h2)12
∂u∂x=x1+x2+h2-a2+c2(x2+h2)2+(a2-c2)-2(a2-c2)(x2-h2)12

2.2编程计算
2.2.1 计算要求
a. 求主剖面上ρs值，（x由-30.0m到+30.0m，步长为1.0m）；
b. h=10,30,50m,ρ1=1.0Ω∙m,ρ2=100.0Ω∙m；
c. c=3.0,5.0,8.0(m) a=0.005,0.01,0.05(m),要求按顺序对应搭配；
d. 将u与∂u∂x的表达式编成二个函数子程序言
e. 以ρs（或ρsρ1）为纵轴，x为横轴，打印出ρs曲线。
2.2.1 C程序源代码
#include
#include
#define r1 1.0 //ρ1
#define r2 100.0 //ρ2
#define h 10 //h=10,30,50m
#define c 3.0 //c= 3.0,5.0,8.0(m)
#define a 0.005 //a=0.005,0.01,0.05(m)
void main()
{
int x;
double rs[61][2]={0}; //ρs
double u[61]={0},px[61]={0};
FILE *fp;
fp=fopen("rs.txt", "w");
double *fun_u(double u1[]);
double *fun_px(double px1[]);
fun_u(u);fun_px(px);
for(x=-30;x<31;x++)
{
rs[x+30][1]=r1*(1+(2*(r1-r2)*a*c)/((a-c)*(r2*a+r1*c))*(1-sqrt((c*c+u[x+30])/(a*a+u[x+30]))
-(1.0/2)*(a*a-c*c)/(sqrt(c*c+u[x+30])*pow(a*a+u[x+30],1.5))*px[x+30]*x));
rs[x+30][1]/=r1;
rs[x+30][0]=x;
fprintf(fp,"%f %f\n",rs[x+30][0],rs[x+30][1]);
printf("%f\n",rs[x+30][1]);
}
fclose(fp);

}
double *fun_u(double u1[])
{
int x;
for(x=-30;x<31;x++)
u1[x+30]=(1.0/2)*(h*h+x*x-a*a-c*c)+(1.0/2)*sqrt(pow(x*x+h*h,2)+pow(a*a-c*c,2)-2*(a*a-c*c)*(x*x-h*h));
return u1;
}
double *fun_px(double px1[])
{
int x;
for(x=-30;x<31;x++)
px1[x+30]=x*(1+(x*x+h*h-a*a+c*c)/sqrt(pow(x*x+h*h,2)+(a*a-c*c)-2*(a*a-c*c)*(x*x-h*h)));
return px1;
}
2.2.1 计算结果
a. 当H=10m,c=3.0m,a=0.005m时，所得结果如图3所示；
b. 当H=30m,c=5.0m,a=0.01m时，所得结果如图4所示；
c. 当H=50m,c=8.0m,a=0.05m时，所得结果如图5所示。

3 结论及几点认识
a. 当ρ2<ρ1，即脉状体为低阻时，ρs或ρs/ρ1在其主剖面的中心上方有极小值。因此，根据ρs曲线主极值点的坐标，可以确定球心在地面的投影位置。
b. 从探测原理上，中间梯度法非常适合直立的高阻异常体。
c. 利用电法勘探寻找矿产资源，是一种简便、快速、经济的方法，经过资料解释和推断，能够大致确定工区的资源分布情况，应综合其他物探资料和当地的地质资料，作出更加准确的解释。
参考文献：
[1]程荃等.地球物理与勘查技术与工程·课程设计指导讲义[M]. 成都理工大学信息工程学院，2005，25～26.
[2] 肖宏跃，雷宛. 地电学教程[M]. 北京：地质出版社，2008，44～46.