复杂形体重磁异常正演.docx

直立长方体模型如图2-1建立三维直角坐标系，x轴指向北方向，y轴指向东方向，z轴铅垂向下。在xy平面下存在一规模为2a*2b*(h2-h1)直立长方体，其在各坐标轴上的投影起始坐标为：1、2- x 坐标的起点和终点1、2- x 坐标的起点和终点1、2- x 坐标的起点和终点，即h1、h2如下图所示：图2- 1 直立长方体 直立长方体的重力异常正演公式直立长方体引起的重力异常表达式为：gr=Gln+r+ln+r-arctanrx1x2y1y2z1z2在上式中，采用CGSM制，重力异常gr单位为Gal，长度单位为cm，剩余密度单位为gcm3，引力常量G为6.67210-8cm3(g.s2)。直立长方体的磁力异常正演公式直立长方体引起的重力异常表达式为：T(x,y,z）=04Mk1lnr+(-x)+ k2lnr+(-y)+ k3lnr+(-z)+ k4arctan-x（-y）（-x）2+r-z+（-z）2+ k5arctan-x（-y）（-y）2+r-z+（-z）2+ k6arctan-x（-y）r-zx1x2y1y2z1z2其中，L0、M0、N0及、分别为地磁场及总磁化强度的方向余弦；则有：k1=M0+N0，k2=L0+N0，k3=L0+M0，k4=L0，k5=M0，k6=-N0。同样，在上式中，采用CGSM制，导磁系数0=1。输入文件格式1.1 输入参数文件格式设计输入参数文件用于记录场的分量标识、导数的阶数、观测面的形态标识、场源参数文件名、测点文件名、输出文件名等信息。本次实验使用的一个参数文件样例cmd.dat为：请选择需要计算的场分量：1 重力或磁力异常 2 重力或磁力异常的导数 1请输入需要计算的导数阶次：1请选择观测面的形态：1 平面规则网 2 平面非规则网3 曲面规则网 4 曲面非规则网1场源参数文件名：source.dat测点位置文件名（规则网时为grd文件，非规则网为xyz文件）：BXYZL.GRD重力异常输出文件：anomaly_gra.GRD磁力异常输出文件（倾斜磁化异常和化极磁异常）：anomaly_mag.GRDanomaly_pol.GRD1.2 场源参数文件格式设计本次实验的场源参数保存在“source.dat ”中。第一列为剩余密度(g/cm3)；第二列为磁化强度(10-6CGSM)；第三列为磁化方向倾角(DEG) ；第四列为磁化方向与x 轴的夹角(DEG) ；第五列为磁化方向与y 轴的夹角(DEG) ；第六列第七列为x 坐标的起点和终点(km) ；第八列第九列为y 坐标的起点和终点(km) ；第十列第十一列为z坐标的起点和终点(km ，向下为正)。如下所示：0.7,20000,50,85,5,-9,-3,4,8,2,70.8,34000,50,85,5,0,7,0,5,3,70.9,17000,50,85,5,-4,8,-8,-3,2,71.3 计算点输入坐标格式设计程序能对如下四种测网进行正演计算：1 平面规则网 2 平面非规则网3 曲面规则网 4 曲面非规则网其中，规则网计算点坐标保存在grd文件中文件后缀为grd；非规则网计算点坐标保存在xyz文件中，文件后缀为txt。对于平面测网，采用只读取第一点高程数据的办法来节省读取外部文件的时间。注：测点坐标文件的后缀必须严格按照如上规定，否则程序将报错并要求重新修改参数文件。 输出文件格式计算完成后程序将自动输出重力异常、磁力异常以及化极磁力异常到外部文件，规则网数据输出到grd文件中，非规则网数据输出到xyz(txt)文件中。PROGRAM Complex_body_forwardIMPLICIT NONECHARACTER*80 cmdfile,sourcefile,stationfile,output_grafile, &output_magfile,output_polfileINTEGER component,order,plane !场的分量标识，导数的阶数，观测面的形态标识INTEGER n_source,n_station,i,j,mpoint,line,cs,cxyzREAL xmin,xmax,ymin,ymaxREAL,ALLOCATABLE:SOURCE(:,:),XYZ(:,:)REAL,ALLOCATABLE:ANOM_GRA(:),ANOM_MAG(:),ANOM_POL(:)cmdfile=cmd.datcs=11cxyz=3CALL INPUT_CMDFILE(component,order,plane,cmdfile,sourcefile, &stationfile,output_grafile,output_magfile,output_polfile)CALL CHECK_CMD(component,order,plane,stationfile)CALL GET_NUMBER_SOURCE(sourcefile,n_source)ALLOCATE(SOURCE(1:n_source,1:cs)CALL INPUT_SOURCE(sourcefile,SOURCE,n_source,cs)CALL GET_n_station(stationfile,plane,xmin,xmax,ymin,ymax, &mpoint,line,n_station)ALLOCATE(XYZ(1:n_station,1:cxyz)CALL INPUT_STATION(stationfile,plane,xmin,xmax,ymin,ymax, &mpoint,line,XYZ,n_station)ALLOCATE(ANOM_GRA(1:n_station),ANOM_MAG(1:n_station)ALLOCATE(ANOM_POL(1:n_station)CALL CAL_ANOM(SOURCE,n_source,cs,XYZ,n_station,cxyz, &ANOM_GRA,ANOM_MAG,ANOM_POL,component)CALL OUTPUT(plane,XYZ,ANOM_GRA,ANOM_MAG,ANOM_POL,n_station,xmin &,xmax,ymin,ymax,mpoint,line,output_grafile,output_magfile &,output_polfile)ENDPROGRAM!读取场的分量标识、导数的阶数、观测面的形态标识!读取场源参数文件名、测点文件名、输出文件名SUBROUTINE INPUT_CMDFILE(component,order,plane,cmdfile,sourcefile, &stationfile,output_grafile,output_magfile,output_polfile)IMPLICIT NONEINTEGER component,order,planeCHARACTER*80 cmdfile,sourcefile,stationfile,output_grafile, &output_magfile,output_polfile OPEN(12,FILE=cmdfile,ACTION=READ)READ(12,*) output_polfileREAD(12,*) output_polfileREAD(12,*) output_polfileREAD(12,*) componentREAD(12,*) output_polfileREAD(12,*) orderREAD(12,*) output_polfileREAD(12,*) output_polfileREAD(12,*) output_polfileREAD(12,*) output_polfileREAD(12,*) output_polfileREAD(12,*) planeREAD(12,*) output_polfileREAD(12,*) sourcefileREAD(12,*) output_polfileREAD(12,*) stationfileREAD(12,*) output_polfileREAD(12,*) output_grafileREAD(12,*) output_polfileREAD(12,*) output_magfileREAD(12,*) output_polfileCLOSE(12)ENDSUBROUTINE!检查输入参数是否正确SUBROUTINE CHECK_CMD(component,order,plane,stationfile)IMPLICIT NONEEXTERNAL UPPERINTEGER component,order,plane,ccCHARACTER*80 stationfile,UPPERCHARACTER signPRINT*,请确认以下信息是否正确：SELECT CASE(component)CASE(1) PRINT*,1 计算重力或磁力异常CASE(2) PRINT*,1 计算重力或磁力异常的导数CASE DEFAULT PRINT*,错误：场分量类型不合法，请修改参数文件！ STOPENDSELECTIF(component=2.AND.order=97.AND.ASC=122) THEN ASC=ASC-32 string(i:i)=CHAR(ASC) ENDIFENDDOUPPER=stringRETURNENDFUNCTION!读取场源的个数SUBROUTINE GET_NUMBER_SOURCE(filename,n)IMPLICIT NONEINTEGER nCHARACTER*80 filenameCHARACTER*20 aaOPEN(101,file=filename,action=read)n=0DO READ(101,*,END=100,ERR=100) aa n=n+1ENDDO100CLOSE(101)ENDSUBROUTINE!读取场源参数SUBROUTINE INPUT_SOURCE(sourcefile,SOURCE,m,n)IMPLICIT NONEINTEGER m,n,i,jREAL ax,ay,inclineCHARACTER*80 sourcefileREAL SOURCE(1:m,1:n)OPEN(12,FILE=sourcefile,ACTION=READ)READ(12,*)(SOURCE(i,j),j=1,n),i=1,m)CLOSE(12)DO i=1,m ax=SOURCE(i,4);ay=SOURCE(i,5);incline=SOURCE(i,3) IF(ax+ay)=90.OR.(ax+ay)=270.OR.ABS(ax-ay)=90.)THEN ELSE PRINT*,错误：磁化方向参数不正确！请修改场源参数文件！ PRINT*,出错行：,i STOP ENDIF IF(incline90)THEN PRINT*,错误：磁倾角参数不正确！请修改场源参数文件！ PRINT*,出错行：,i STOP ENDIFENDDOENDSUBROUTINE!读取测量点的个数SUBROUTINE GET_n_station(stationfile,plane,xmin,xmax,ymin,ymax, &mpoint,line,n_station)IMPLICIT NONEINTEGER plane,mpoint,line,n_stationREAL xmin,xmax,ymin,ymaxCHARACTER*80 stationfileIF(plane=1.OR.plane=3) THEN CALL GET_GRID(stationfile,mpoint,line,xmin,xmax,ymin,ymax) n_station=mpoint*lineELSE mpoint=0;line=0;xmin=0;xmax=0;ymin=0;ymax=0 CALL GET_NUMBER_SOURCE(stationfile,n_station)ENDIFENDSUBROUTINE!读取规则网测网网格参数SUBROUTINE GET_GRID(inputfile,mpoint,line,xmin,xmax,ymin,ymax)IMPLICIT NONECHARACTER*80 inputfileINTEGER mpoint,lineREAL xmin,xmax,ymin,ymaxOPEN(101,FILE=inputfile,ACTION=READ)READ(101,*)READ(101,*)mpoint,line,xmin,xmax,ymin,ymaxCLOSE(101)ENDSUBROUTINE!读取测点坐标SUBROUTINE INPUT_STATION(stationfile,plane,xmin,xmax,ymin,ymax, &mpoint,line,XYZ,n_station)IMPLICIT NONEREAL xmin,xmax,ymin,ymax,dx,dy,tempINTEGER plane,mpoint,line,n_stationREAL XYZ(1:n_station,1:3)CHARACTER*80 stationfileREAL i,jOPEN(12,FILE=stationfile,ACTION=READ)IF(plane=1.OR.plane=3) THEN DO i=1,5,1 READ(12,*) ENDDO dx=(xmax-xmin)/(mpoint-1) dy=(ymax-ymin)/(line-1) DO j=1,line DO i=1,mpoint XYZ(j-1)*mpoint+i,1)=(i-1)*dx+xmin XYZ(j-1)*mpoint+i,2)=(j-1)*dy+ymin ENDDO ENDDO IF(plane=1) THEN READ(12,*) XYZ(1,3) XYZ(:,3)=XYZ(1,3) ELSE READ(12,*) XYZ(:,3) ENDIFELSEIF(plane=2.OR.plane=4) THEN IF(plane=2) THEN READ(12,*) XYZ(1,:) XYZ(:,3)=XYZ(1,3) DO i=2,n_station READ(12,*) XYZ(i,1),XYZ(i,2) ENDDO ELSE DO i=1,n_station READ(12,*) XYZ(i,:) ENDDO ENDIFELSEPRINT*,测网类型不合法！ENDIFENDSUBROUTINE!计算各测点的异常值SUBROUTINE CAL_ANOM(SOURCE,n_source,cs,XYZ,n_station,cxyz, &ANOM_GRA,ANOM_MAG,ANOM_POL,component)IMPLICIT NONEEXTERNAL Dg,ElementDg,DTREAL Dg,ElementDg,dg0,DT,dt0INTEGER n_source,cs,n_station,cxyz,i,j,componentREAL SOURCE(1:n_source,1:cs),XYZ(1:n_station,1:cxyz)REAL ANOM_GRA(1:n_station),ANOM_MAG(1:n_station)REAL ANOM_POL(1:n_station)ANOM_GRA=0.;ANOM_MAG=0.;ANOM_POL=0.IF(component=1) THEN DO i=1,n_station,1 DO j=1,n_source dg0=Dg(XYZ(i,1)*1E5,XYZ(i,2)*1E5,XYZ(i,3)*1E5,SOURCE(j,6:7)*1E5 & ,SOURCE(j,8:9)*1E5,SOURCE(j,10:11)*1E5,SOURCE(j,1), & ElementDg) ANOM_GRA(i)=ANOM_GRA(i)+dg0 dt0=DT(XYZ(i,1)*1E5,XYZ(i,2)*1E5,XYZ(i,3)*1E5,SOURCE(j,6:7)*1E5 & ,SOURCE(j,8:9)*1E5,SOURCE(j,10:11)*1E5,SOURCE(j,2)*1E-6, & SOURCE(j,3),SOURCE(j,4),SOURCE(j,5) ANOM_MAG(i)=ANOM_MAG(i)+dt0 dt0=DT(XYZ(i,1)*1E5,XYZ(i,2)*1E5,XYZ(i,3)*1E5,SOURCE(j,6:7)*1E5 & ,SOURCE(j,8:9)*1E5,SOURCE(j,10:11)*1E5,SOURCE(j,2)*1E-6, & 90.,SOURCE(j,4),SOURCE(j,5) ANOM_POL(i)=ANOM_POL(i)+dt0 ENDDO ENDDOELSEIF(component=2) THEN PRINT*,重磁异常导数计算程序正在开发中. PRINT*,$目前暂不支持此项功能$ STOPENDIFENDSUBROUTINE!*计算某源点在计算点引起的重力异常程序集开始*FUNCTION Dg(X,Y,Z,X_Source,Y_Source,Z_Source,Density,ElementDg)IMPLICIT NONEREAL X,Y,Z,X_Source(1:2),Y_Source(1:2),Z_Source(1:2), &Density,ElementDgREAL Dg,ELMTINTEGER I,J,K,SIGNDg=0ELMT=0DO I=1,2DO J=1,2DO K=1,2SIGN=(-1)*(I+J+K)ELMT=ElementDg(X,Y,Z,X_Source(I), &Y_Source(J),Z_Source(K),Density)Dg=Dg+ELMT*SIGNENDDOENDDOENDDOENDFUNCTION!计算积分表达式中的某一分量FUNCTION ElementDg(X,Y,Z,X_Source,Y_Source,Z_Source,Density)IMPLICIT NONEEXTERNAL ATANXY,LOG0REAL G,ATANXY,LOG0PARAMETER(G=6.672E-8)REAL ElementDg,X,Y,Z,X_Source,Y_Source,Z_Source,Density,R,T1,T2,T3R=SQRT(X_Source-X)*2+(Y_Source-Y)*2+(Z_Source-Z)*2)T1=(X_Source-X)*LOG0(Y_Source-Y)+R)T2=(Y_Source-Y)*LOG0(X_Source-X)+R)T3=(Z_Source-Z)*ATANXY(X_Source-X)*(Y_Source-Y),(Z_Source-Z)*R)ElementDg=-G*Density*(T1+T2-T3)ENDFUNCTION!*计算某源点在计算点引起的重力异常程序集结束*!*计算某源点在计算点引起的磁力异常程序集开始*FUNCTION DT(x,y,z,X_Source,Y_Source,Z_Source,moment,incline,ax,ay)IMPLICIT NONEEXTERNAL ATANXY,LOG0REAL x,y,z,moment,incline,ax,ay,DT,ATANXY,LOG0REAL X_Source(1:2),Y_Source(1:2),Z_Source(1:2)INTEGER i,j,k,signREAL hax,hay,za,mx,my,mz,pi,r,xs,ys,zspi=3x=moment*COSD(incline)*COSD(ax)my=moment*COSD(incline)*COSD(ay)mz=moment*SIND(incline)hax=0.;hay=0.;za=0.DO i=1,2,1 DO j=1,2,1 DO k=1,2,1 sign=(-1)*(i+j+k) xs=X_Source(i);ys=Y_Source(j);zs=Z_Source(k) r=SQRT(x-xs)*2+(y-ys)*2+(z-zs)*2) hax=hax+sign*(-mx*ATANXY(xs-x)*(ys-y),(xs-x)*2+r*(zs-z)+ & (zs-z)*2)+my*LOG0(r+zs-z)+mz*LOG0(r+ys-y) hay=hay+sign*(mx*LOG0(r+zs-z)-my*ATANXY(xs-x)*(ys-y),( & (ys-y)*2+r*(zs-z)+(zs-z)*2)+mz*LOG0(r+xs-x) za=za+sign*(mx*LOG0(r+ys-y)+my*LOG0(r+xs-x)-mz* & ATANXY(xs-x)*(ys-y),r*(zs-z) ENDDO ENDDOENDDODT=hax*COSD(incline)*COSD(ax)+hay*COSD(incline)*COSD(ay)+ & za*SIND(incline)DT=DT/4/piENDFUNCTION!*计算某源点在计算点引起的磁力异常程序集结束*!计算atan(x/y)FUNCTION ATANXY(x,y)IMPLICIT NONEREAL x,y,ATANXY,eps,pi2pi2=1.5707963267948966192313216916398eps=4./HUGE(x)IF(ABS(x)=eps.AND.ABS(y)eps) THEN