1、一、基本原理 2二、输入、输出格式设计 .21、输入数据： 22、输出数据： 2三、总体设计 3四、测试结果 41、观测面重力异常图： 42、计算面重力异常图： 4五、结论及建议 5附录：源程序代码.5一、基本原理本次实验的基本原理为根据位场特征选择一个离散的、等效的解析表达式，并推导该解析表达式的各种处理和转换表达式，利用已知条件建立线性方程组，并求解该线性方程组得到等效源的有关参数。这里采用分块迭代的方案进行。把得到的等效源参数代入前面的等效源解析式就可以进行各种处理和转换重力异常计算公式：= (-)(-)2+(-)2+(-)232二、输入、输出格式设计1、输入数据：观测面 x 坐标（向东
2、）为 X1y 坐标（向北）为 Y1z 坐标（铅垂向下）为 Z1重力异常为 G1等效源 x 坐标（向东）为 X3，y 坐标（向北）为 Y3，z 坐标（铅垂向下）为Z3，重力异常为 G3等效源与观测面点的距离为 R1，与计算面点的距离为 R2,密度为 P2、输出数据：计算面 x 坐标（向东）为 X2y 坐标（向北）为 Y2z 坐标（铅垂向下）为 Z2重力异常为 G2三、总体设计程序流程图：开始输入有关参数观测面坐标值及重力值，计算面坐标值从文件中读取已知点个数、计算点个数布置等效源，建立方程组运用 LU 分解，求解方程，得到等效源参数计算计算面重力异常结束四、测试结果1、观测面重力异常图：-261
3、840-6182-26840-6182- -40-图 4.1 观测面重力异常图2、计算面重力异常图：-261840-6182-61840-2618- -40-五、结论及建议为了给平面异常转换处理和反演方法的应用创造条件，起伏地形上重力异常转换的最重要的内容，是由起伏地形上重力异常换算出某一平面上的重力异常，简称为“曲化平” 。其中，等效源的建立及方程求解是关键。编制程序时，要注意循环，谨慎处理，以防出错。附录：源程序代码!观测面 x 坐标（向东）为 X1，y 坐标（向北）为 Y1，z 坐标（铅垂向下）为 Z1，重力异常为 G1!计算面 x 坐标（向东）为 X2，y 坐标（向北）为 Y2，z 坐
4、标（铅垂向下）为 Z2，重力异常为 G2!等效源 x 坐标（向东）为 X3，y 坐标（向北）为 Y3，z 坐标（铅垂向下）为 Z3，重力异常为 G3!等效源与观测面点的距离为 R1，与计算面点的距离为 R2,密度为 P图 4.2 计算面重力异常图program dengxiaoyuanparameter(G=0.0667)integer n1,n2real,allocatable:X1(:),X2(:),X3(:),Y1(:),Y2(:),Y3(:),Z1(:),Z2(:),&Z3(:),G1(:),G2(:),G3(:),P(:),P1(:),R1(:,:),R2(:,:),A(:,:),L
5、(:,:),&U(:,:)real dz,c,d,e,fcall duqugeshu(gravity.dat,n1)call duqugeshu(xyz.dat,n2)allocate(X1(1:n1),X2(1:n2),X3(1:n1),Y1(1:n1),Y2(1:n2),Y3(1:n1),&Z1(1:n1),Z2(1:n2),Z3(1:n1),G1(1:n1),G2(1:n2),G3(1:n1),P(1:n1),&R1(1:n1,1:n1),R2(1:n2,1:n2),A(1:n1,1:n1),L(1:n1,1:n1),U(1:n1,1:n1)&,P1(1:n1)call duqushuj
6、u(gravity.dat,xyz.dat,X1,Y1,Z1,G1,X2,Y2,Z2,n1,&n2)call dengyuan(G1,n1,P,R1,X1,X3,Y1,Y3,Z1,Z3,G)do i=1,n1do j=1,n1A(i,j)=G*(Z1(i)-Z3(j)/(R1(i,j)*3)end doend docall lufenjie(A,L,U,P,R1,G1,G,Z1,Z3,n1)call jisuan(R2,X2,Y2,Z2,G2,X3,Y3,Z3,G,P,n1,n2)call output(X2,Y2,G2,n2,out.dat)deallocate(X1,X2,X3,Y1,Y2
7、,Y3,Z1,Z2,Z3,G1,G2,G3,P,R1,R2,A,L,U,P1)end!读入个数subroutine duqugeshu(filename,number)integer numbercharacter*(*) filenameOpen(11,file=filename,status=old)Number=0Do while(.not.eof(11)Read (11,*,end=100,ERR=100)x,y,zNumber=number+1100 End doclose(11)end subroutine!读入数据subroutine duqushuju(filename1,fi
8、lename2,X1,Y1,Z1,G1,X2,Y2,Z2,n1,&n2)character*(*) filename1,filename2integer n1,n2real X1(1:n1),Y1(1:n1),Z1(1:n1),G1(1:n1),X2(1:n2),Y2(1:n2),&Z2(1:n2)open(22,file=filename1)do i=1,n1read(22,*) X1(i),Y1(i),Z1(i),G1(i)end doopen(33,file=filename2)do i=1,n2read(33,*) X2(i),Y2(i),Z2(i)end doend subrouti
9、ne!等效源subroutine dengyuan(G1,n1,P,R1,X1,X3,Y1,Y3,Z1,Z3,G)integer n1real Greal G1(1:n1),P(1:n1),R1(1:n1,1:n1),X1(1:n1),X3(1:n1),Y1(1:n1),&Y3(1:n1),Z1(1:n1),Z3(1:n1),G11(1:n1)do i=1,n1X3(i)=X1(i)Y3(i)=Y1(i)Z3(i)=-10.0end dodo i=1,n1do j=1,n1R1(i,j)=sqrt(X1(i)-X3(j)*2+(Y1(i)-Y3(j)*2+&(Z1(i)-Z3(j)*2)end
10、 doend dodo i=1,n1 do j=1,n1G11(i)=G11(i)+G*P(j)*(Z1(i)-Z3(j)/(R1(i,j)*3)end do end dodo i=1,n1G1(i)=G11(i)end doend subroutine!LU 分解subroutine lufenjie(A,L,U,P,R1,G1,G,Z1,Z3,n1)real G1(1:n1),P(1:n1),R1(1:n1,1:n1),A(1:n1,1:n1),L(1:n1,1:n1),&U(1:n1,1:n1),B(1:n1,1:n1),P1(1:n1),Z1(1:n1),Z3(1:n1)integer
11、 n1real dz,G,c,d,e,fdo i=1,n1do j=1,n1if(i=j) thenL(i,j)=1.0end ifif(ji) thenL(i,j)=0.0end ifif(ij) thenU(i,j)=0.0end ifend doend dodo j=1,n1U(1,j)=A(1,j)end dodo i=2,n1L(i,1)=A(i,1)/U(1,1)end dodo r=2,n1do i=r,n1do k=1,r-1c=c+L(r,k)*U(k,i)end doU(r,i)=A(r,i)-cend doend dodo r=2,n1do i=r+1,n1if(r/=n
12、1) thendo k=1,r-1d=d+L(i,k)*U(k,r)end doL(i,r)=(A(i,r)-d)/U(r,r)end ifend doend doP1(1)=G1(1)do i=2,n1do k=1,i-1e=e+L(i,k)*P1(k)end doP1(i)=G1(i)-eend doP(n1)=P1(n1)/U(n1,n1)do i=n1-1,1,-1do k=i+1,n1f=f+U(i,k)*P(k)end doP(i)=(P1(i)-f)/U(i,i)end doend subroutine!计算平面重力值subroutine jisuan(R2,X2,Y2,Z2,G
13、2,X3,Y3,Z3,G,P,n1,n2)integer n1,n2real Greal X2(1:n2),X3(1:n1),Y2(1:n2),Y3(1:n1),Z2(1:n2),Z3(1:n1),&G2(1:n2),P(1:n1),R2(1:n2,1:n2),G22(1:n2)do i=1,n2do j=1,n1R2(i,j)=sqrt(X2(i)-X3(j)*2+(Y2(i)-Y3(j)*2+&(Z2(i)-Z3(j)*2)end doend dodo i=1,n2 do j=1,n1G22(i)=G22(i)+G*P(j)*(Z2(i)-Z3(j)/(R2(i,j)*3)end do end dodo i=1,n2G2(i)=G22(i)end doend subroutine!输出subroutine output(X2,Y2,G2,n2,filename)character*(*) filenameinteger n2real X2(1:n2),Y2(1:n2),G2(1:n2)open(77,file=filename)do i=1,n2write(77,*) X2(i),Y2(i),G2(i)end doclose(77)end subroutine
