物探数据处理实验报告

本科生实验报告实验课程物探数据处理学院名称地球物理学院专业名称勘查技术与工程学生姓名00000学生学号*************指导教师李勇实验地点5417实验成绩二〇一七年九月二〇一七年十月物探数据处理实验报告实验一1.1实验目的本次实验的目的是将课本上给出的地质体情况，和计算参量的代数公式用c 语言通过程序表现出来。

在得知已知参量，例如形状大小、空间位置的情况下将其带入编写好的程序，从而得出地下板状体的物性值数据，并将所得数据其通过成图数据成图后在用反演原理与提前已知的异常体情况对比，判断所得结果的正误，并分析最终的实验结果。

1.2实验原理下图1-1为实验原理示意图。

图1-1 二维板模型利用如下公式进行c语言程序的实现：1.3实验参数中点横坐标x0(m)：1000 中点深度坐标z0(m)：1000 板高度(m)：400 板长度(m)：200 磁化强度：2000 二维板角度：45、90 磁化角度：90 测点k：100}1.4 实验结果图示图1-2 磁异常X分量图1-3 磁异常Z分量图1-4 重力异常1.5 实验结果分析根据实验数据结果图可知，通过观察x方向的磁分量可以看出：图中存在着明显的两个相互对称的正负异常，可知在地质体的上部积累了一定的负磁荷，而在地质体的下部积累了同样大小的正磁荷，实验中磁化角度取90度，二维板倾角取不同的角度。

通过观察z方向上的磁异常可知：地质体呈现出中部正异常的对称图形，总体呈现先增后减的趋势。

通过观察重力异常的图像可知：得出的结果基本与已知相符，所以通过重力异常、磁场x方向异常、磁场z方向异常数据曲线，就可对所求数据做相应的分析。

实验二2.1 实验目的通过给定的已知相关物性参数，建立截面为多边形的水平柱体模型，带入所编写好的c 语言程序中得到对应的正演结果。

主要考查了c 语言程序进编程，和对所得结果的理解，分析，判断能力。

2.2 实验原理图2-1为实验原理图图2-1 多边形截面水平柱体截面形状根据下列公式进行编程得到实验数据，作图分析。

；]1111[)1(2)1(2)11)(1(221212ln )1(2)1(2)11)(1(21{12x i zi tg x i z i tg x i x i z i z i z i x i z i x i x i x i z i x i z i x i x i x i z i z i z i x i z i x i z i z i n i G g ++---⋅-++-++-+-++++++⋅-++-++-+-+⋅=∑=∆σ；)(21Q M z P M x X +=∆π ；)(21P M z Q M x Z -=∆π式中：；z ix i z i x i x i x i z i z i z i z i x i z i tg x i z i tg x i x i z i z i x i x i z i z i n i Q 221212ln )1(2)1(2)1(221]1111[)1(2)1(2)1)(1({1++++⋅+-+-+-+⋅-++---⋅+-+-++--+=∑=；z ix i z i x i x i x i z i z i x i x i z i z i x i z i tg x i z i tg x i x i z i z i z i z i n i P 221212ln )1(2)1(2)1()1(21]1111[)1(2)1(2)1(2{1++++⋅+-+-++--+⋅+++---⋅+-+-+-+=∑=2.3 实验参数多边形半径r:200m ；磁化强度M:2000 拟合边数:10 磁化角度ct:90,120 测点k:100 中点深度z0(m):1000 中点横坐标x0(m):10002.3.1 讨论abtg 1-的值ab实际为由坐标原点对多边形第i 边的夹角，在具体计算中应当进行讨论： 当0.115->a 时，如果a b >0.0；当a>0.0,b>0.0,则θ=a btg 1-；当a<0.0，b<0.0，则θ=—π+abtg 1-；如果a b <0.0，则θ=π+ab tg 1-；当0.115-≤a 时，若b<0.0，则θ=—2π；若b>=0,则θ=2π2.4 实验结果图2-2 重力异常图2-3 磁异常X分量图2-4 磁异常Z分量2.5 实验结果分析通过重力异常图可知，通过观察x方向的磁分量可以看出：图中仍然存在着明显的两个相互对称的正负异常，故我们仍然可以认为地质体的上部积累了一定的负磁荷，而在地质体的下部积累了同样大小的正磁荷。

所以左端出现了较大的正异常，右端出现了较大的正异常。

通过观察z方向上的磁异常可知：地质体在中线（x=50）呈现出中部正异常的对称图形，总体呈现先增后减的趋势。

通过观察重力异常的图像可知：得出的结果基本与已知相符，出现了一个重力异常值曲线，最大值点对应的测线坐标与地质体实际位置一致。

实验三3.1 实验目的本次实验的目的是通过给定的地下异常体参量，对地下的正长方体进行正演计算，然后将计算参量的代数公式用c 语言通过程序表现出来。

在得知已知参量，例如形状大小、空间位置的情况下将其带入编写好的程序，并通过积分来实现正演。

从而得出地下板状体的物性值数据，并将所得数据其通过成图数据成图后在用反演原理与提前已知的异常体情况对比，判断所得结果的正误，并分析最终的实验结果。

将其运用的主要方法有通过划分单元来进行体积分，即8个体积元叠加即为立方体所对应的积分值，偶数为正，奇数为负，则可以计算出相应重磁异常。

3.2 实验原理利用如下公式进行编程：；2020202020201ln ln ({c z c z b y b y a x ax }zk)R(z yk)xk)(y (x tg zk)(z R]xk)[(x yk)(y R]yk)[(y xk)x G Δg -+-+-+----⋅-++--++---=σ；2020********yk)]}-(y ln[R M z zk)]-(z ln[R M y x k)R-(x zk)-yk)(z -(y tg 1-M x {-41c z cz b y b y a x a x Δx -+-+-+++++=π；2020********x k)]}-(x ln[R M z ]yk)R-(y zk)-x k)(z -(x tg 1-M y -zk)]-(z ln[R M x {41c z c z b y b y a x a x Δy -+-+-+++++=π；2020********}zk)R-(z yk)-x k)(y -(x tg 1-M z x k)]-(x ln[R M y -yk)]-(y ln[R M x {41c z c z b y b y a x a x Δz -+-+-+-++=π式中：])(2)(2)(2[2/1zk z yk y xk x R -+-+-=；3.3 实验参数列数:101 行数: 101 磁化强度:2000 磁倾角:45 磁偏角:45 立方体上坐标:（300,300,10） 立方体下坐标:（700,700,410） 密度:1 起终坐标:（0,1000）3.4 实验结果图示图3-1 重力异常平面图 图3-2 磁异常x 分量平面图3-3 磁异常y 分量平面 图3-4 磁异常z 分量平面图3.5实验结果分析根据实验数据结果图可知，通过观察x方向的磁分量可以看出：图中存在着两个对称的正负异常，且异常呈现出相似性形态幅值基本一致。

通过观察z方向上的磁异常可知：地质体呈现出中部正异常的对称图形，切呈现出对角线对称的形态上方为正下方为负的负异常。

通过观察重力异常的图像可知：得出的结果基本与已知相符，得到一个明显的类似同心圆的重力异常，边界较为清晰。

实验程序1#include<stdio.h>#include<math.h>#define PI 3.1415926#define G 6.67e-6void main(){FILE *fp1,*fp2,*fp3;fp1=fopen("g.txt","w");fp2=fopen("deltax.txt","w");fp3=fopen("deltaz.txt","w");double r1,r2,r3,r4; //定义板截面角点A,B,C,D到P点之间的距离double f1,f2,f3,f4; //定义r1.r2.r3.r4与X轴正向夹角int k,i; //定义P点的坐标xk,zk;double xk,zk;double cgm=2.67e+3;double x0=1000.0; //中点横坐标double z0=1000.0;double l=400.0; //板高度double b=200.0; //板长度double I,M=2000.0;double z1,z2,g,deltax,deltaz;double J=90.0*PI/180.0; //计算二维板角度I=90.0*PI/180.0; //计算磁化角度double a[100]={0.0};double b1[100]={0.0};double ct[100]={0.0};for(k=0;k<101;k++){xk=k*20.0;zk=0.0;b1[0]=z0-zk-l*sin(J);a[0]=xk-x0+b+l*cos(J);b1[1]=z0-zk+l*sin(J);a[1]=xk-x0+b-l*cos(J);b1[2]=z0-zk-l*sin(J);a[2]=xk-x0-b+l*cos(J);b1[3]=z0-zk+l*sin(J);a[3]=xk-x0-b-l*cos(J);r1=sqrt(pow(b1[0],2.0)+pow(a[0],2.0));r2=sqrt(pow(b1[1],2.0)+pow(a[1],2.0));r3=sqrt(pow(b1[2],2.0)+pow(a[2],2.0));r4=sqrt(pow(b1[3],2.0)+pow(a[3],2.0));for(i=0;i<4;i++){if(fabs(a[i])>1.0e-15){if((b1[i]/a[i])>0.0){if(a[i]>0.0&&b1[i]>0.0){ct[i]=atan(b1[i]/a[i]);}else if(a[i]<0.0&&b1[i]<0.0){ct[i]=atan(b1[i]/a[i])-PI;}}else if((b1[i]/a[i])<0.0){ct[i]=PI+atan(b1[i]/a[i]);}}else if(fabs(a[i])<=1.0e-15){if(b1[i]<0.0){ct[i]=-PI/2.0;}else{ct[i]=PI/2.0;}}f1=PI-ct[0];f2=PI-ct[1];f3=PI-ct[2];f4=PI-ct[3];z1=z0-l*sin(J);z2=z1+l*sin(J);g=2.0*G*cgm*((z2*(f2-f4)-z1*(f1-f3))+xk*(sin(J)*sin(J)*log((r2*r3)/(r 1*r4))+cos(J)*sin(J)*(f1-f2-f3+f4))+2.0*b*(sin(J)*sin(J)*log(r4/r3)+c os(J)*sin(J)*(f3-f4)));deltax=((M*sin(J))/(2.0*PI))*(log((r2*r3)/(r1*r4))*cos(J-I)-(sin(J-I) )*(f1-f2-f3+f4));deltaz=((M*sin(J))/(2.0*PI))*((sin(J-I))*log((r2*r3)/(r1*r4))+(cos(J-I))*(f1-f2-f3+f4));fprintf(fp1,"%f\n",g); //将所得数据写入文件fprintf(fp2,"%f\n",deltax);fprintf(fp3,"%f\n",deltaz);}}fclose(fp1);fclose(fp2);fclose(fp3);实验程序2include"stdio.h"#include"math.h"#define PI 3.1415926 //圆周率#define G 6.67e-11 // 万有引力常数#define r 200.0 //多边形半径#define M 2000.0 //磁化强度#define thegama 2.67e+3void main(){double ct11(double a,double b); //声明子函数FILE *fp1,*fp2,*fp3;fp1=fopen("deltag.txt","w");fp2=fopen("deltax1.txt","w");fp3=fopen("deltaz1.txt","w");double x0=1000.0;double z0=1000.0;int n=10,N=100;double R=200.0;double M=2000.0;double cgm=2.67e+3;double ct=PI/4.0;double Mx=M*cos(ct);double Mz=M*sin(ct);double x[10000]={0.0},x1[1000]={0.0};double z[10000]={0.0},z1[1000]={0.0};double xk[10000]={0.0},zk[1000]={0.0};double sum1=0.0;double Q=0.0;double P=0.0;doublea1,a2,a3,a4,a5,a6,a7,aa,bb,ct1,ct2,deltag,deltax1,deltaz1,a,b;int i,k;for(i=0;i<n+1;i++)/////////////////////////∠坐标{x1[i]=x0+R*cos(i*2*PI/n);z1[i]=z0-R*sin(i*2*PI/n);}for(k=0;k<N;k++){xk[k]=20.0*k;///////////////////平面坐标zk[k]=0.0;}for(k=0;k<N;k++){sum1=0.0;Q=0.0;P=0.0;for(i=0;i<n+1;i++){x[i]=xk[k]-x1[i];z[i]=zk[k]-z1[i];}for(i=0;i<n;i++){a1=(z[i+1]-z[i])*(x[i]*z[i+1]-x[i+1]*z[i])/(pow(z[i+1]-z[i],2.0)+ pow(x[i+1]-x[i],2.0));a2=(x[i+1]-x[i])*(x[i]*z[i+1]-x[i+1]*z[i])/(pow(z[i+1]-z[i],2.0)+pow( x[i+1]-x[i],2.0));a3=(z[i+1]-z[i])*(x[i]-x[i+1])/(pow(z[i+1]-z[i],2.0)+pow(x[i]-x[i +1],2.0));a4=pow(z[i+1]-z[i],2.0)/(pow(z[i+1]-z[i],2.0)+pow(x[i]-x[i+1],2.0));aa=pow(x[i+1],2.0)+pow(z[i+1],2.0);bb=pow(x[i],2.0)+pow(z[i],2.0);ct1=ct11(z[i],x[i]); //调用子函数ct2=ct11(z[i+1],x[i+1]);a5=(a1/2.0)*log(aa/bb)+a2*(atan(ct1)-atan(ct2));sum1+=a5;a6=a3*(atan(ct1)-atan(ct2))-(a4/2.0)*log(aa/bb);Q+=a6;a7=a4*(atan(ct1)-atan(ct2))+(a3/2.0)*log(aa/bb);P+=a7;}deltag=2.0*G*cgm*sum1;deltax1=(Mx*P+Mz*Q)/(2.0*PI);deltaz1=(Mx*Q-Mz*P)/(2.0*PI);fprintf(fp1,"%lf\n",deltag);fprintf(fp2,"%lf\n",deltax1);fprintf(fp3,"%lf\n",deltaz1);}fclose(fp1);fclose(fp2);fclose(fp3);}///////////////////////////////////////////double ct11(double a,double b) //判断子函数{double ct;if(fabs(a)>1.0e-15){if((b/a)>0.0){if(a>0.0&&b>0.0)ct=atan(b/a);elsect=atan(b/a)-PI;}if((b/a)<0.0){ct=PI+atan(b/a);}}else{if(b<0.0){ct=-PI/2.0;}else{ct=PI/2.0;}}return(ct);}实验程序3#include <stdlib.h>#include <stdio.h>#include <math.h>#define l 101 // 列数#define h 101 // 行数#define M 2000.0 //磁化强度#define II 60 // 磁倾角#define aa 60 //磁偏角#define PI 3.1415926 //圆周率#define G 6.674E-11 // 万有引力常数#define SI 1e6 // 将m/s^2换算为g.u.的比例因子#define DENSITY 1E3 // 将g/cm^3换算为kg/m^3的比例因子//计算atan值double batan(double m,double n,double x,double r);//Kernel函数,用途：为Cuboid函数准备double Kernel(double dXk, double dYk, double dZk,//dXk, dYk, dZk: 观测点坐标double dXs, double dYs, double dZs);//dYs, dZs, dZs: 立方体坐标，s取1,2 //Cuboid函数,用途: 计算立方体重力异常,参数说明double Cuboid(double dXk, double dYk, double dZk,//dXk, dYk, dZk: 观测点坐标double dX1, double dY1, double dZ1,// dX1, dY1, dZ1: 立方体坐标1double dX2, double dY2, double dZ2,//dX2, dY2, dZ2: 立方体坐标2double dDens); //dDens: 立方体密度// Cix函数double Cix(double dXk, double dYk, double dZk,double dXs, double dYs, double dZs,double Mx, double My, double Mz);// Cixxx函数: 用于计算立方体△x异常double Cixxx(double dXk, double dYk, double dZk,double dX1, double dY1, double dZ1,double dX2, double dY2, double dZ2,double Mx, double My, double Mz );// Ciy函数double Ciy(double dXk, double dYk, double dZk,double dXs, double dYs, double dZs,double Mx, double My, double Mz);// Ciyyy函数: 用于计算立方体△Y磁异常double Ciyyy(double dXk, double dYk, double dZk,double dX1, double dY1, double dZ1,double dX2, double dY2, double dZ2,double Mx, double My, double Mz );// Ciz函数double Ciz(double dXk, double dYk, double dZk,double dXs, double dYs, double dZs,double Mx, double My, double Mz);// Cizzz函数: 用于计算立方体△Z磁异常double Cizzz(double dXk, double dYk, double dZk,double dX1, double dY1, double dZ1,double dX2, double dY2, double dZ2,double Mx, double My, double Mz );void main(){FILE *f1,*f2,*f3,*f4; //定义文件指针f1 = fopen("重力异常.txt", "w"); //打开文件"重力异常.txt"f2 = fopen("磁异常X分量.txt", "w"); //打开文件"磁异常X分量.txt"f3 = fopen("磁异常Y分量.txt", "w"); //打开文件"磁异常Y分量.txt" f4 = fopen("磁异常Z分量.txt", "w"); //打开文件"磁异常Z分量.txt"int i, j;// 网格参数double dMinX = 0.0, dMaxX = 1000.0; // x坐标范围double dMinY = 0.0, dMaxY = 1000.0; // y坐标范围double dDX = (dMaxX - dMinX) / (l-1); // x方向上的网格间距double dDY = (dMaxY - dMinY) / (h-1); // y方向上的网格间距double px[101]={0.0},py[101]={0.0}; // x,y坐标构成的向量double pz[101][101]; // z坐标构成的矩阵double dg[101][101]={0.0}; // 定义重力异常数组double dx[101][101]={0.0},dy[101][101]={0.0},dz[101][101]={0.0}; //定义磁异常X分量、磁异常Y分量以及磁异常Z分量数组double I=II*PI/180; //角度转换double a=aa*PI/180; //角度转换// 指定坐标范围for (i = 0; i < h; i++){py[i] = dMinY + i*dDY;for (j = 0; j < l; j++){pz[i][j] = 0.0;}}for(i=0; i<l; i++){px[i] = dMinX + i*dDX;}// 定义立方体参数double dX1 = 300, dY1 = 300, dZ1 = 10; // 坐标1double dX2 = 700, dY2 = 700, dZ2 = 410; // 坐标2double dDens = 1.0;// 密度double Mx = M*cos(I)*cos(a);double My = M*cos(I)*sin(a);double Mz = M*sin(I); //磁分量// 生成重力与磁分量异常for (i = 0; i < h; i++){for (j = 0; j < l; j++){dg[i][j] = Cuboid(px[i], py[j], pz[i][j],dX1, dY1, dZ1,dX2, dY2, dZ2, dDens);dx[i][j] = Cixxx(px[i], py[j], pz[i][j],dX1, dY1, dZ1,dX2, dY2, dZ2, Mx, My, Mz);dy[i][j] = Ciyyy(px[i], py[j], pz[i][j],dX1, dY1, dZ1,dX2, dY2, dZ2, Mx, My, Mz);dz[i][j] = Cizzz(px[i], py[j], pz[i][j],dX1, dY1, dZ1,dX2, dY2, dZ2, Mx, My, Mz);fprintf(f1,"%d\t%d\t%f\n",i*10,j*10,dg[i][j]);//将所得到的重力异常数据写入文件fprintf(f2,"%d\t%d\t%f\n",i*10,j*10,dx[i][j]);//将所得到的磁异常X分量数据写入文件 fprintf(f3,"%d\t%d\t%f\n",i*10,j*10,dy[i][j]);//将所得到的磁异常Y分量数据写入文件 fprintf(f4,"%d\t%d\t%f\n",i*10,j*10,dz[i][j]);//将所得到的磁异常Z分量数据写入文件}}}//arctan的子程序double batan(double m,double n,double x,double r){double q;if(fabs(x*r)>1.0e-15){if((m*n/(x*r))>0.0){if((m*n)>0.0&&(x*r)>0.0) q=atan(m*n/(x*r));else q=-PI+atan(m*n/(x*r));}else q=PI+atan(m*n/(x*r));}else{if((m*n)<0.0) q=-PI/2;else q=PI/2;}return q;}// Kernel函数double Kernel(double dXk, double dYk, double dZk,double dXs, double dYs, double dZs){double dX=dXs-dXk, dY=dYs-dYk, dZ=dZs-dZk;double dR=sqrt(dX*dX+dY*dY+dZ*dZ);if(dZ==0){return -dX*log(dR+dY)-dY*log(dR+dX);}else{return -dX*log(dR+dY)-dY*log(dR+dX)+dZ*atan(dX*dY/(dZ*dR));}}// Cuboid函数: 用于计算立方体重力异常double Cuboid(double dXk, double dYk, double dZk,double dX1, double dY1, double dZ1,double dX2, double dY2, double dZ2,double dDens){return SI*DENSITY*G*dDens*(-Kernel(dXk, dYk , dZk, dX1, dY1, dZ1) //SI取数据为毫伽，见常数定义+Kernel(dXk, dYk , dZk, dX1, dY1, dZ2) //8个体积叠加即为立方体所对应的积分值+Kernel(dXk, dYk , dZk, dX1, dY2, dZ1) //偶数为正，奇数为负-Kernel(dXk, dYk , dZk, dX1, dY2, dZ2)+Kernel(dXk, dYk , dZk, dX2, dY1, dZ1)-Kernel(dXk, dYk , dZk, dX2, dY1, dZ2)-Kernel(dXk, dYk , dZk, dX2, dY2, dZ1)+Kernel(dXk, dYk , dZk, dX2, dY2, dZ2));}// Cix函数double Cix(double dXk, double dYk, double dZk,double dXs, double dYs, double dZs,double Mx, double My, double Mz){double dX=dXs-dXk, dY=dYs-dYk, dZ=dZs-dZk;double dR=sqrt(dX*dX+dY*dY+dZ*dZ);return My*log(dR+dZ)+Mz*log(dR+dY)-Mx*batan(dY,dZ,dX,dR);}// Cixxx函数: 用于计算立方体△X磁异常double Cixxx(double dXk, double dYk, double dZk,double dX1, double dY1, double dZ1,double dX2, double dY2, double dZ2,double Mx, double My, double Mz ){return 1.0/(4*PI)*(-Cix(dXk, dYk , dZk, dX1, dY1, dZ1, Mx, My,Mz) +Cix(dXk, dYk , dZk, dX1, dY1, dZ2, Mx, My,Mz)+Cix(dXk, dYk , dZk, dX1, dY2, dZ1, Mx, My,Mz) -Cix(dXk, dYk , dZk, dX1, dY2, dZ2, Mx, My,Mz)+Cix(dXk, dYk , dZk, dX2, dY1, dZ1, Mx, My,Mz)-Cix(dXk, dYk , dZk, dX2, dY1, dZ2, Mx, My,Mz)-Cix(dXk, dYk , dZk, dX2, dY2, dZ1, Mx, My,Mz)+Cix(dXk, dYk , dZk, dX2, dY2, dZ2, Mx, My,Mz));}// Ciy函数double Ciy(double dXk, double dYk, double dZk,double dXs, double dYs, double dZs,double Mx, double My, double Mz){double dX=dXs-dXk, dY=dYs-dYk, dZ=dZs-dZk;double dR=sqrt(dX*dX+dY*dY+dZ*dZ);return Mx*log(dR+dZ)+Mz*log(dR+dX)-My*batan(dX,dZ,dY,dR);}// Ciyyy函数: 用于计算立方体△Y磁异常double Ciyyy(double dXk, double dYk, double dZk,double dX1, double dY1, double dZ1,double dX2, double dY2, double dZ2,double Mx, double My, double Mz ){return 1.0/(4*PI)*(-Ciy(dXk, dYk , dZk, dX1, dY1, dZ1, Mx, My,Mz) +Ciy(dXk, dYk , dZk, dX1, dY1, dZ2, Mx, My,Mz)+Ciy(dXk, dYk , dZk, dX1, dY2, dZ1, Mx, My,Mz) -Ciy(dXk, dYk , dZk, dX1, dY2, dZ2, Mx, My,Mz)+Ciy(dXk, dYk , dZk, dX2, dY1, dZ1, Mx, My,Mz)-Ciy(dXk, dYk , dZk, dX2, dY1, dZ2, Mx, My,Mz)-Ciy(dXk, dYk , dZk, dX2, dY2, dZ1, Mx, My,Mz)+Ciy(dXk, dYk , dZk, dX2, dY2, dZ2, Mx, My,Mz));}// Ciz函数double Ciz(double dXk, double dYk, double dZk,double dXs, double dYs, double dZs,double Mx, double My, double Mz){double dX=dXs-dXk, dY=dYs-dYk, dZ=dZs-dZk;double dR=sqrt(dX*dX+dY*dY+dZ*dZ);return Mx*log(dR+dY)-My*log(dR+dX)-Mz*batan(dX,dY,dZ,dR);}// Cizzz函数: 用于计算立方体△Z磁异常double Cizzz(double dXk, double dYk, double dZk,double dX1, double dY1, double dZ1,double dX2, double dY2, double dZ2,double Mx, double My, double Mz ){return 1.0/(4*PI)*(-Ciz(dXk, dYk , dZk, dX1, dY1, dZ1, Mx, My,Mz) +Ciz(dXk, dYk , dZk, dX1, dY1, dZ2, Mx, My,Mz)+Ciz(dXk, dYk , dZk, dX1, dY2, dZ1, Mx, My,Mz) -Ciz(dXk, dYk , dZk, dX1, dY2, dZ2, Mx, My,Mz)+Ciz(dXk, dYk , dZk, dX2, dY1, dZ1, Mx, My,Mz)-Ciz(dXk, dYk , dZk, dX2, dY1, dZ2, Mx, My,Mz)-Ciz(dXk, dYk , dZk, dX2, dY2, dZ1, Mx, My,Mz)+Ciz(dXk, dYk , dZk, dX2, dY2, dZ2, Mx, My,Mz}。