《复杂形体重磁异常正演.docx》由会员分享，可在线阅读，更多相关《复杂形体重磁异常正演.docx（20页珍藏版）》请在金锄头文库上搜索。
1、直立长方体模型如图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)。直立长
2、方体的磁力异常正演公式直立长方体引起的重力异常表达式为：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 输入参数文件格式设计输入参数文件用于记
3、录场的分量标识、导数的阶数、观测面的形态标识、场源参数文件名、测点文件名、输出文件名等信息。本次实验使用的一个参数文件样例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.
4、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,
5、-3,2,71.3 计算点输入坐标格式设计程序能对如下四种测网进行正演计算：1 平面规则网 2 平面非规则网3 曲面规则网 4 曲面非规则网其中，规则网计算点坐标保存在grd文件中文件后缀为grd；非规则网计算点坐标保存在xyz文件中，文件后缀为txt。对于平面测网，采用只读取第一点高程数据的办法来节省读取外部文件的时间。注：测点坐标文件的后缀必须严格按照如上规定，否则程序将报错并要求重新修改参数文件。 输出文件格式计算完成后程序将自动输出重力异常、磁力异常以及化极磁力异常到外部文件，规则网数据输出到grd文件中，非规则网数据输出到xyz(txt)文件中。PROGRAM Complex_bod
6、y_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
7、(:),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
8、 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)ALLOCA
9、TE(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!读取场的分量标识、导数的阶数、观测面的形态标识!读取场源参数文件名、
10、测点文件名、输出文件名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=REA
11、D)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_pol
12、fileREAD(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 NONEEX
13、TERNAL 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.order0) THEN PRINT*,错误：导数阶次不合法，请修改参数文件！ STO
14、PENDIFSELECT CASE(plane)CASE(1) PRINT*,2 平面规则网 CASE(2) PRINT*,2 平面非规则网CASE(3) PRINT*,2 曲面规则网CASE(4) PRINT*,2 曲面非规则网CASE DEFAULT PRINT*,错误：测网类型不合法，请修改参数文件！ STOPENDSELECTcc=LEN_TRIM(stationfile)IF(plane=1.OR.plane=3).AND. &(UPPER(stationfile(cc-3:cc).ne.GRD) THEN PRINT*,错误：测点文件类型不合法，请修改参数文件！ STOPENDIFIF(plane=2.OR.plane=4).AND. &(UPPER(stationfile(cc-3:cc).ne.TXT) THE
