实验三多次迭代滑动平均法程序设计姓 名：专业：勘查技术与工程学 号：指导教师：纪新林.基本原理:1.多次迭代法的计算公式:G(0)(i, j) = G (i, j)艺 艺 G(p)(i + k, i +1)G(P+1) (i, j) = +=Nrk = Mr艺艺1l=-Nr l=-Mr式中：G(i, j)(i = l,m; j = l,n)为总场值;G(0)(i, j)(i = 1,m; j = 1,n)为迭代的初始值；P为迭代次数，Nr为窗口沿y轴的半径，Mr为窗口沿x轴的半径。2.多次迭代法可解决的问题：可以进行场的分离，得到深部区域场和剩余场，同时也可以消除随机异常。二输入输出数据格式设计：X方向窗口半径：MradiusY方向窗口半径：Nradius 迭代次数：Interationmax 输入总场文件名：Totalfilename 输出区域场文件名：Regfilename 输出剩余场文件名：Resfilename 总场值：G 区域场值：GReg 剩余场值：GRes 点数：mpoint线数：nline点坐标的最小值：Xmin点坐标的最大值：Xmax线坐标的最小值：Ymin线坐标的最大值：Ymax三总体设计I：输入有关参数：文件名，迭代次数，窗口半径 输入点线数，X，Y坐标的最大值和最小值 输入总场P：循环进行滑动平均计算0：输出区域场和剩余场四测试结果1. 测试参数如下:(1) 网格数据保存在Totalfilename.grd中。网格数据按grd格式保存，有27 条线、27个点，点(线)坐标从-26m26m，点线距为2m。(2) 迭代次数3,窗口半径为6m,相关参数保存在CMD. txt文件中：Tot alfilename.grdRegfilename.grdResfilename.grd3662. 测试结果：得到区域场和剩余场的两个网格化文件：Regf il ename.grd ，Resf il ename.grd 用Surfer 8分别绘制等值线图：分别得到总场等值线图，区域场等值线图，剩余场等值线图：2520151050-5 -10 -15 -20 -25-25-20-15-10-50510152025总场等值线图0.38 .36 0.34 0.32 0.3 0.28 0.26 0.24 0.220.20.180.160.140.120.10.080.060.040.0200.10.0950.090.0850.080.0750.070.0650.060.0550.050.0450.040.0350.03-25-20-15-10-50510152025区域场等值线图2520151050-5-10-15-20-25-25-20-15-100.280.260.240.220.20.180.160.140.120.10.080.060.040.02 0 -0.02-0.04-0.06-50510152025剩余场等值线五结论及建议:多次迭代滑动平均法可以进行场的分离，将一个重力异常分离为深部的区域 异常和剩余场，同时能够消除随机虚假异常。一般来说，窗口越大，滑动平均值反映的区域异常越深；迭代次数越多，得 到的区域场越平滑，得到的剩余场越明显，反复迭代基本上可以消除随机异常和 假异常，达到压制干扰的效果。本次测试迭代了三次，由以上图可知，迭代得到的区域异常比总异常平滑了 很多，而剩余场基本上与总场一致，有细微的差别。总体来说得到了实验要求的 效果，基本上反映了多次迭代滑动平均法进行场的分离的规律。为了使结果更加明显，可以增加迭代次数。附录：源程序代码程序说明mulfi_iter_mean -重力异常数据的多次迭代滑动平均处理输入参数:Mradius - X方向窗口半径Nradius - Y方向窗口半径Interationmax -迭代次数mpoint -点数nline -线数Xmin - X测线最小值所用文件:所用文件:中间变量:中间变量:输出参数:输出参数:Totalfilename -输入总场文件名 Xmax - X测线最大值 Regfilename -输出区域场文件名 Ymin - Y测线最小值 Resfilename -输出剩余场文件名 Ymax - Y测线最大值 G -总场值CMD.txt -输入相关参数Totalfilename.grd -总场输入文件Regfilename.grd -区域场输出文件Resfilename.grd -剩余场输出文件number -中间整型变量GReg -区域场值sum -中间实型变量GRes -剩余场值Gmax -网格点处异常数据的最大值Gmax -网格点处异常数据的最小值Regfilename.grd -区域场输出文件Resfilename.grd -剩余场输出文件:input_CMD-输入相关文件名和数据 read_mpointnline -从输入文件中读入点线个数和X，Y测线最大值最小值inputGRD -输入总场值到数组G mulfi_iter_mea -滑动平均法 mean_sub -迭代过程 outputGRD -输出网格化文件作者：韩自强!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc!滑动平均法的主程序PROGRAM mulfi_iter_meanCHARACTER*80 Totalfilename,Regfilename,ResfilenameINTEGER mrradius,nradiusINTEGER interationmaxINTEGER mpoint,nlineREAL xmin,xmax,ymin,ymaxREAL,ALLOCATABLE: G(:,:),GReg(:,:),GRes(:,:)CALL input_CMD(Totalfilename,Regfilename,Resfilename,interationmax,mradius,nradius)CALL read_mpointnline(Totalfilename,mpoint,nline,xmin,xmax,ymin,ymax)ALLOCATE(G(l:mpoint,l:nline),GReg(l:mpoint,l:nline),GRes(1:mpoint,1:nline)CALL inputGRD(Totalfilename,mpoint,nline,G)CALL mulfi_iter_mea(mpoint,nline,Gmradius,nradius,interationmax,GReg,GRes)CALL outputGRD(Regfilename,mpoint,nline,GReg,xmin,xmax,ymin,ymax)CALL outputGRD(Resfilename,mpoint,nline,GRes,xmin,xmax,ymin,ymax) DEALLOCATE(GGReg,GRes)*END PROGRAM mulfi_iter_meanSUBROUTINE input_CMD(Totalfilename,Regfilename,Resfilename,interationmax,mradius,nradius)CHARACTER*80 Totalfilename,Regfilename,ResfilenameINTEGER interationmax,mrradius,nradiusOPEN(20,file=CMD.txt)READ(20,*) TotalfilenameREAD(20,*) RegfilenameREAD(20,*) ResfilenameREAD(20,*) interationmaxREAD(20,*) mradiusREAD(20,*) nradiusclose(20)END SUBROUTINE input_CMD!从输入文件中读入点线个数和X, Y测线最大值最小值SUBROUTINE read_mpointnline(Totalfilename,mpoint,nline, xmin,xmax,ymin,ymax)CHARACTER*80 TotalfilenameINTEGER mpoint,nlineREAL xmin,xmax,ymin,ymaxOPEN(20,file=Totalfilename)READ(20,*)READ(20,*) mpoint,nlineREAD(20,*) xmin,xmaxREAD(20,*) ymin,ymaxCLOSE(20)END SUBROUTINE read_mpointnline!输入总场值到数组GSUBROUTINE inputGRD(Totalfilename,mpoint,nline,G)CHARACTER*80 TotalfilenameINTEGER mpoint,nlineREAL G(l:mpoint,l:nline)OPEN(20,file=Totalfilename)READ(20,*)READ(20,*) mpoint,nlineREAD(20,*)READ(20,*)READ(20,*)READ(20,*) (G(I,J),I=l,mpoint),J=l,nline)CLOSE(20)END SUBROUTINE inputGRD!滑动平均法SUBROUTINE mulfi_iter_mea(mpoint,nline,G,mradius,nradius,*interationmax,GReg,GRes)INTEGER mrradius,nradiusINTEGER interationmaxINTEGER mpoint,nlineREAL G(l:mpoint,l:nline),GReg(l:mpoint,l:nline),* GRes(l:mpoint,l:nline)GRes=GDO K=l,interationmaxCALL mean_sub(mpoint,nline,mradius,nradius,GReg,GRes)GRes=GRegEND DOGRes=G-GRegEND SUBROUTINE mulfi_iter_mea!迭代过程SUBROUTINE mean_sub(mpoint,nli