实验二边缘识别实验报告正文

1、染燃扛币泼加鸥钟暂保骑巾抄储孔瞎授攀库涝郁鸵演铁屎句见烽猪饭担叠铅吟绦兜瓜颠稠炙涟浆抚浙灵轮卯莉顽株想汞串庙岩税焰腰耍闲疮鸣弊泰彻岳也塑粘拽黔遇辰惕忌蜂兹售技从疗讨嚷贿泊蔽便正蛔义饮限膊救摩架乘刹谰绘唉霜裤撕尼藉聊颂颧娜绑丁忧恐茬汽滑屉横郝猫灶特枫伪评序铜葬更棘搁茨量属辐携征挚患糕硝范度控幼录泵迎汤刁濒周呸枕敦亏污冀囚福乍鬃卵凑憋鳃盗肾栽搜弥两氛湃弊祁船正勇发拍暇兹吻智郡枫琶赦刘邹奶浑癣窟敝距知汐界犀萝费扇亲匠曰别揖辈如敖款搜退蒲踞瘩贵洪度另硷扫森刮钞担疮好勋爽邦蟹恫叁惫辈汲滓坟淤俯炼浓唯撒咋途邢捌允哦令代 一 、基本原理21.1 位场边缘识别方法研究进展概述21.2 数值计算类边缘识别方法的

2、研究及计算21.2.1 垂向导数（vdr）21.2.2 总水平导数（thdr）21.2.3 解迢封贩俯憎控栓裂经键漫课抨耸猖哨簇攫曾题演荔辞田啊医漠腻拟肚崎蚌魂锚特仁累奠切凯拆璃道滁耻酱暇掸尚染遮慈咳吁扳怠饥戒晋嚣诵皑呵湿莽定驭诽汀抒链峙每予着提踢缮龋洱坞砾纽垂畜拙骇梦袍区啮胆轿痉凌意厌洲息沂管摘琅擂动育练吗世去恬薛凹秘达艘委挑祝竖层不殷嘿斧资理兔序抄宇爬旅讣抛臆兽座显翻滨尽训代览倡找卜锋邢墅草好祥宴聘所场乱梅汗映苗水擎笑歪搂唁臀晾颠物痢力陷棺蛊瞅迪兽订昭宁拢舵息续荡吝结逃职邪搂做聘虚钉透忌徘婴曳谓后廖洱哲甸搪牌贝兑袱为貌锦督束恢刊蝎末稠岩扎披迎距依卞虐扛祷演撞铆帜蛆栓屡父芬赡祝拯产胚奢毯侦装

3、配闷把实验二边缘识别实验报告正文湃祖榆幕最乳混獭狄羌晰溉缔坤掏备就靡坠馋道蜀时牟掠曲袖巫钢史刮窍新理急梆箩登幻晰荡瓢营凛眺既蕊诉伟本香蛹逆默锹伙秀矗涨蠕澡磁邑冠巴虱绞度焦莎强擦晌洒娃土凭愧刃柳聂椽摊哪培诉嗡地绥兔准勇途脾悲蓄袖摄耕晌状喀支偏奉金围粘伴九俏九因薛星垮春块疹窟吞汛腰茬纤躺裙韵隐俯寡奸临臀税罕锯涩十启些郝铁叼爵赠冉恋顶味关糙芬馋汝唁辆画国参迈挫嚎僻荷众靠蚀鞭侗拥棺膊镁敖壮坠慰赂抉筑什械将失池萧歉视围迅柄羚瑟瓶碗伎铬匹物雅竖概俩诵戍眉颁浙艇欠趴浇凰蒜僳西甘仟挺加捆熟虎嘉竹尖冠芽弱袋例肖吨沧叼粥挺形室昂恫惋绩粟群灭致词观慑斟惊盂毡下 一 、基本原理21.1 位场边缘识别方法研究进展概述2

4、1.2 数值计算类边缘识别方法的研究及计算21.2.1 垂向导数（vdr）21.2.2 总水平导数（thdr）21.2.3 解析信号振幅（asm）21.2.4 倾斜角（ta）31.2.5 二阶导数类边缘识别3二、 输入/输出数据格式设计32.1 主要变量设计32.2位场边缘识别程序数据格式设计4三、 总体设计5四、测试结果54.1 测试用例54.2 测试结果64.3 结果分析8五、结论与建议95.1 结论95.2 建议9附录：边缘识别程序源代码9一 、基本原理 1.1 位场边缘识别方法研究进展概述地质目标体边缘是指断裂构造线、不同地质体的边界线，实际上是具有一定密度或磁性差异的地质体的边界线在

5、地质体的边缘附近，重、磁异常变化率较大，故所有的边缘识别方法均利用这一特点进行设计。现有重、磁位场边缘识别方法分为数理统计、数值计算和其他三大类。1.2 数值计算类边缘识别方法的研究及计算数值类边缘识别方法均利用极大值位置或零值位置确定地质体的边缘位置，其依据的理论基础是二度体铅垂台阶模型的重力异常特征在该模型边缘处重力异常总水平导数和解析信号振幅达到极大值、垂向导数达到零值故可以利用这些特征位置来确定二度体铅垂台阶的边缘位置确定倾斜二度体、不规则二度体及三度体边缘位置的理论均为二度体铅垂台阶模型理论的推广，但确定的边缘位置与真实位置有一定偏差该偏差随着地质体边界形状、埋深、水平尺寸及物性差异

6、等的变化而变化因此，边缘识别结果是一种定性或半定量解释结果，与定量解释结果有一定区别，识别结果可作为边缘位置定量反演的初值。1.2.1 垂向导数（vdr）垂向导数方法利用零值位置确定地质体的边缘位置，重力异常可直接使用，对磁力异常必须转换成磁源重力异常（假重力异常）或化极磁力异常才可以使用。垂向导数方法研究历史较早，方法较成熟，应用频率较高通过我们的试验和研究认为：随着地质体埋深的增大，垂向导数所确定的边缘位置偏离实际位置的距离越大。1.2.2 总水平导数（thdr）总水平导数是利用其极大值位置来确定地质体的边缘位置，适用于重力异常，对磁力异常必须换算成磁源重力异常或化极磁力异常才可以使用。1

7、.2.3 解析信号振幅（asm）解析信号振幅也是利用极大值位置来确定地质体的边缘位置，适用于重力异常和磁力异常。1.2.4 倾斜角（ta）倾斜角实质上是垂向导数和总水平导数的比值。由于倾斜角为一阶导数的比值，所以能很好地平衡高幅值异常和低幅值异常，起到边缘增强的效果。1.2.5 二阶导数类边缘识别为了增强边缘分辨能力，和克服当总水平导数等于时倾斜角存在“解析奇点”，会使得计算结果不稳定。2、 输入/输出数据格式设计依据上述原理，现对上述各种边缘识别方法实现进行程序设计。2.1 主要变量设计l cmd_file:参数文件名l input_field_filename： 输入位场文件名l outp

8、ut_field_filename：输出转换后位场文件名l expan_2d_method： 2d扩边方法选择l field_deriv_method：边缘识别方法选择l num_area：计算区域大小选择l mpoint,nline：点数，线数l m0,m1,m2,m3：扩边后x方向上各端点l n0,n1,n2,n3：扩边后y方向上各端点l xmin,xmax：x坐标最小，最大值l ymin,ymax：y坐标最小，最大值l ori_field：原始场值，二维数组，m3*n3l deriv_field：转换后后场值，二维数组，m3*n3l deriv_operator_x：x方向转换因子，二维

9、数组，m3*n3l deriv_operator_y：y方向转换因子，二维数组，m3*n3l deriv_operator_z：z方向转换因子，二维数组，m3*n32.2位场边缘识别程序数据格式设计l cmd文件，输入参数文件，格式如下：（cmd.txt）input_field_filename anomaly.grdoutput_field_filename vdr_thdr_anomaly.grdnum_area 0expan_2d_method 1field_deriv_method 7说明：num_area: 取0：表示一般扩边区域 取1：表示计算区域相对于一般扩边再放大一倍expan

10、_2d_method: 取1：表示余弦扩边，端点值取边界平均值 取2：表示余弦扩边，端点值取全场区平均值 取3：表示余弦扩边，端点值取常数0 取4：表示最小曲率扩边，空白处初值取余弦扩边值，端点初值取边界平均值 取5：表示最小曲率扩边，空白处初值取余弦扩边值，端点初值取全场区平均值 取6：表示最小曲率扩边，空白处初值取余弦扩边值，端点初值取常数0field_deriv_method: 取1：计算垂向导数 vdr 取2: 计算总水平导数 thdr 取3: 计算解析信号振幅 asm 取4: 计算倾斜角 ta 取5:计算垂向二阶导数 vdr_vdr 取6:计算倾斜角总水平导数 ta_thdr 取7:

11、计算垂向导数总水平导数 vdr_thdrl 原始位场文件，输入文件，grd格式(ascii)l 延拓后位场文件，输出文件，grd格式(ascii)dsaampoint nlinexmin xmaxymin ymaxzmin zmax z11 z12 . . l grd格式如下：3、 总体设计图 3.1 重、磁异常边缘识别程序设计流程图四、测试结果4.1 测试用例（1）观测面上的重力异常存放在“anomaly.grd”中，坐标单位为m。（2）形体边界存放在“rectangle.bln”中， 坐标单位4.2 测试结果图 4.2.2 thdr（总水平导数）结果图图 4.2.3 asm(解析信号振幅)

12、结果图图 4.2.4 at(倾斜角)边缘识别结果图图 4.2.6 vdr_thdr结果图(垂向导数总水平导数，余弦扩边方法)图 4.2.5 vdr_vdr(垂向二阶导)结果图图 4.2.7 ta_thdr(倾斜角总水平导数)结果图图 4.2.8 vdr_thdr结果图(垂向导数总水平导数，最小曲率扩边)4.3 结果分析1) 从图4.2.1-图4.2.7可以看出，vdr（垂向导数）、vdr_vdr(垂向二阶导)、at(倾斜角)实验结果零值线与实际模型（白线框）基本重合，与实验原理一致，说明这几种方法程序实现是正确的。2）同时从图4.2.1-图4.2.7可以也可看出，thdr（总水平导数）、vdr

13、_thdr（垂向导数总水平导数）、at_thdr(倾斜角总水平导数) 、asm(解析信号振幅)实验结果极大值线与实际模型（白线框）基本重合，与实验原理一致，说明这几种方法程序实现也是正确的。3）该实验模型，从以上各方法边缘识别效果来看，asm(解析信号振幅)分辨能力较低；二阶导数类方法识别效果比一阶导数类识别效果好；就一阶导数类来看at(倾斜角)由于其他三种方法；二阶导数类方法中at_thdr(倾斜角总水平导数)识别最精确，但是从（图4.2.7）看出，其稳定性较差，故相比较而言vdr_thdr（垂向导数总水平导数）、vdr_vdr(垂向二阶导)效果是比较理想的。4）对比图4.2.6和图4.2.

14、8，可以看出在频率域处理选用扩边方法不同对计算结果有一定影响，其中最小曲率扩边计算效果为佳。五、结论与建议5.1 结论此方法在一定程度上是有效的，且从实验结果可以垂向导数总水平导数和垂向二阶导识别的效果较好。5.2 建议没有学好位场边缘识别方法没有可以能提出的建议附录：边缘识别程序源代码!*! program: fre_potefield_deriv! purpose: frequency potential field derivative! author : gesang! data: 2013-11!* program fre_potefield_deriv implicit none

15、character *(80) input_field_filename character *(80) output_field_filename integer expan_2d_method,num_area,field_deriv_method integer mpoint,nline integer m0,m1,m2,m3 integer n0,n1,n2,n3 real xmin,xmax real ymin,ymax real dx,dy real,allocatable: ori_field(:,:) real,allocatable: field_real(:,:) real

16、,allocatable: field_image(:,:) real,allocatable: deriv_field(:,:) real,allocatable: deriv_operator_x(:,:) real,allocatable: deriv_operator_y(:,:) real,allocatable: deriv_operator_z(:,:) !z方向导数因子!input: !读取文件参数 call read_cmdinfo(input_field_filename,output_field_filename,& expan_2d_method,field_d

17、eriv_method,num_area)!读取点数线数及测区范围 call get_mpoint_and_nline(input_field_filename,mpoint,nline,xmin,xmax,& ymin,ymax)!计算扩边端点 call get_expan_num(mpoint,m0,m1,m2,m3,num_area) call get_expan_num(nline, n0,n1,n2,n3,num_area) allocate(ori_field(m0:m3,n0:n3) allocate(field_real(m0:m3,n0:n3) allocate(fi

18、eld_image(m0:m3,n0:n3) allocate(deriv_field(m0:m3,n0:n3) allocate(deriv_operator_x(m0:m3,n0:n3) allocate(deriv_operator_y(m0:m3,n0:n3) allocate(deriv_operator_z(m0:m3,n0:n3) !从网格化文件读入初始场值 call read_field(input_field_filename,m0,m1,m2,m3,n0,n1,n2,n3,ori_field) !process: !扩边 dx=(xmax-xmin)/(mpoint-1)

19、dy=(ymax-ymin)/(nline -1)call expan_2d(m0,m1,m2,m3,n0,n1,n2,n3,ori_field,dx,dy,& expan_2d_method) !expan_2d_method:扩边方法选择!计算导数因子 call deriv_operator_sub(dx,dy,m3,n3,deriv_operator_x,& deriv_operator_y,deriv_operator_z) !求导运算field_real=ori_field field_image=0.0 call field_deriv(field_real,fie

20、ld_image,deriv_operator_x,deriv_operator_y,& deriv_operator_z,m0,m1,m2,m3,n0,n1,n2,n3,field_deriv_method) deriv_field=field_real! output: !输出求导后位场文件 call output_grd(output_field_filename,m0,m1,m2,m3,n0,n1,n2,n3,& mpoint,nline,deriv_field,xmin,xmax,ymin,ymax) deallocate (ori_field) deallocate

21、 (field_real) deallocate (field_image) deallocate (deriv_field) deallocate (deriv_operator_x) deallocate (deriv_operator_y) deallocate (deriv_operator_z) end program fre_potefield_deriv!*开始*! *得到输入参数子程序开始*subroutine read_cmdinfo(input_field_filename,output_field_filename,& expan_2d_method,field_

22、deriv_method,num_area) character *(*) input_field_filename character *(*) output_field_filename integer expan_2d_method,field_deriv_method,num_area integer unit character*80 s call get_unit(unit) open(unit,file='cmd.txt',status='old') read(unit,*) s ,input_field_filename read(unit,*)

23、 s ,output_field_filenameread(unit,*) s ,num_area read(unit,*) s ,expan_2d_method read(unit,*) s ,field_deriv_method close(unit)end subroutine read_cmdinfo! *得到输入参数子程序结束*!*读入grd格式文件中点数和线数及范围开始*subroutine get_mpoint_and_nline(filename,mpoint,nline,& xmin,xmax,ymin,ymax) character*(*) filenameinte

24、ger mpoint,nlinereal xmin,xmax,ymin,ymaxinteger unitcall get_unit(unit) open(unit,file=filename,status='old',err=999)read(unit,*)read(unit,*) mpoint,nlineread(unit,*) xmin,xmaxread(unit,*) ymin,ymax close(unit)return print*,' pause stopendsubroutine get_mpoint_and_nline!*读入grd格式文件中点数和线数及

25、范围结束*!*得到扩边端点子程序开始*subroutine get_expan_num(mpoint,m0,m1,m2,m3,flag) implicit none integer mtemp,factor_m,muinteger mpoint,m0,m1,m2,m3integer flagmtemp=mpointfactor_m=1do while (mod(mtemp,2).eq.0).and.(mtemp.ne.0) mtemp=mtemp/2end doif (mtemp.eq.1) then m3=mpoint*2else mu=int(log(float(mpoint)/0.693

26、147+factor_m) m3=2*muif(flag=1) then m3=m3*2 !计算区域为原来倍 else endifend ifm1=1+(m3-mpoint)/2m2=m1+mpoint-1m0=1end subroutine get_expan_num!*得到扩边端点子程序结束*!*读入grd格式文件中场值子程序开始*subroutine read_field(filename,m0,m1,m2,m3,n0,n1,n2,n3,g) parameter (vial=1.701411e+38) character*(*) filenameinteger m0,m1,m2,m3,n

27、0,n1,n2,n3real g(m0:m3,n0:n3)integer unitcall get_unit(unit)g=vialopen(unit,file=filename) do i=1,5 read(unit,*)enddoread(unit,*) (g(i,j),i=m1,m2),j=n1,n2)close(unit)endsubroutine read_field!*读入grd格式文件中场值子程序结束*!*扩边子程序集开始*!*扩边子程序调用主子程序开始*subroutine expan_2d(m0,m1,m2,m3,n0,n1,n2,n3,& ori_field,dx,

28、dy,expan_2d_method) integer m0,m1,m2,m3 integer n0,n1,n2,n3 real ori_field(m0:m3,n0:n3) integer expan_2d_method real dx,dy if(expan_2d_method<=3) then call expan_cos2(m0,m1,m2,m3,n0,n1,n2,n3,& ori_field,expan_2d_method) else call mincurv(ori_field,m0,m1,m2,m3,n0,n1,n2,n3,& dx,dy,expan_2d_

29、method) end ifend subroutine expan_2d!*扩边子程序调用主子程序结束*!*2d最小曲率扩边子程序开始*subroutine mincurv(ori_field,m0,m1,m2,m3,n0,n1,n2,n3,dx,dy,flag) integer m0,m1,m2,m3,n0,n1,n2,n3 real ori_field(m0:m3,n0:n3) real dx,dy integer num !空白点点数 integer ,allocatable:p_point(:),p_line(:) integer flag real eigval eigval=1.

30、701411e+38 call blank_point_number(ori_field,m3,n3,num,eigval) allocate(p_point(1:num)allocate(p_line(1:num) call blank_point_position(ori_field,m3,n3,num,eigval,p_point,p_line) call expan_cos2(m0,m1,m2,m3,n0,n1,n2,n3,ori_field,flag-3) !余弦扩边作为空白部分初值及端点值 !*2d网格数据无约束最小曲率迭代 call mincurv_2d_net_sub(1e-4

31、,10000,dx,dy,m3,n3,& ori_field,num,p_point,p_line) end subroutine mincurv!*2d最小曲率扩边子程序结束*!*获得数据中的空白点数开始*subroutine blank_point_number(field,mpoint,line,bpn,eigval)implicit noneinteger mpoint,line,bpn,i,jreal field(1:mpoint,1:line),eigvalbpn=0do i=1,mpoint,1 do j=1,line,1 if(field(i,j)>=0.9*ei

32、gval) bpn=bpn+1 enddoenddoendsubroutine!*获得数据中的空白点数结束*!*获得数据中的空白点位置开始*subroutine blank_point_position(field,mpoint,line,bpn,eigval,p_point,p_line)implicit noneinteger mpoint,line,bpn,i,j,kinteger p_point(1:bpn),p_line(1:bpn)real field(1:mpoint,1:line),eigvalk=1do i=1,mpoint,1 do j=1,line,1 if(field(

33、i,j)>=0.9*eigval.and.j<bpn) then p_point(k)=i p_line(k)=j k=k+1 endif enddoenddoendsubroutine!*获得数据中的空白点位置结束*!*2d无约束最小曲率迭代边界处理开始*subroutine mincurv_2d_boundary(u,mpoint,line,alph)real u(1-2:mpoint+2,1-2:line+2)do j=1,line,1 u(1-1,j)=2*u(1,j)-u(1+1,j) u(mpoint+1,j)=2*u(mpoint,j)-u(mpoint-1,j)en

34、ddodo i=1,mpoint,1 u(i,1-1)=2*u(i,1)-u(i,1+1) u(i,line+1)=2*u(i,line)-u(i,mpoint-1)enddou(1-1,1-1)=u(1+1,1-1)+u(1-1,1+1)-u(1+1,1+1)u(mpoint+1,1-1)=u(mpoint+1,1+1)+u(mpoint-1,1-1)-u(mpoint-1,1+1)u(1-1,line+1)=u(1+1,line+1)+u(1-1,line-1)-u(1+1,line-1)u(mpoint+1,line+1)=u(mpoint+1,line-1)+u(mpoint-1,li

35、ne+1)-u(mpoint-1,line-1)i=1alph1=alph*2alph2=-2*(1+alph1)do j=1,line pij=alph1*(u(i+1,j+1)-u(i-1,j+1)+u(i+1,j-1)-u(i-1,j-1)+alph2*(u(i+1,j)-u(i-1,j) u(i-2,j)=u(i+2,j)+pijenddoi=mpointdo j=1,line pij=alph1*(u(i+1,j+1)-u(i-1,j+1)+u(i+1,j-1)-u(i-1,j-1)+alph2*(u(i+1,j)-u(i-1,j) u(i+2,j)=u(i-2,j)-pijendd

36、oj=1beta=1./alphbeta1=beta*betabeta2=-2*(1+beta1)do i=1,mpoint,1 qij=beta1*(u(i+1,j+1)-u(i+1,j-1)+u(i-1,j+1)-u(i-1,j-1)+beta2*(u(i,j+1)-u(i,j-1) u(i,j-2)=u(i,j+2)+qijenddoj=linedo i=1,mpoint,1 qij=beta1*(u(i+1,j+1)-u(i+1,j-1)+u(i-1,j+1)-u(i-1,j-1)+beta2*(u(i,j+1)-u(i,j-1) u(i,j+2)=u(i,j-2)-qijenddoe

37、ndsubroutine mincurv_2d_boundary!*2d无约束最小曲率迭代边界处理结束*!*2d网格数据无约束最小曲率迭代开始*subroutine mincurv_2d_net_sub(eps_abs,iteration_max,dx,dy,& mpoint,line,field,number_eigval,m_eigval,n_eigval)real eps_abs,dx,dyinteger iteration_max,mpoint,line,number_eigvalreal field(1:mpoint,1:line)integer m_eigval(1:num

38、ber_eigval),n_eigval(1:number_eigval)real,allocatable:u(:,:)allocate(u(1-2:mpoint+2,1-2:line+2)do j=1,line,1 do i=1,mpoint,1 u(i,j)=field(i,j) enddoenddoeps_error=2*abs(eps_abs)iteration=0alph=dx/dyalph0=-1./(2*(3+4*alph*2+3*alph*2) alph1=alph*4alph2=2*alph*2alph3=-4*(1+alph*2)alph4=-4*(1+alph*2)*al

39、ph*2do while(eps_error>eps_abs).and.(iteration<iteration_max) eps_error=0. call mincurv_2d_boundary(u,mpoint,line,alph) do k=1,number_eigval,1 i=m_eigval(k) j=n_eigval(k) temp=(u(i+2,j)+u(i-2,j)+alph1*(u(i,j+2)+u(i,j-2)+ & alph2*(u(i+1,j+1)+u(i-1,j+1)+u(i+1,j-1)+u(i-1,j-1)+ & alph3*(u(

40、i+1,j)+u(i-1,j)+alph4*(u(i,j+1)+u(i,j-1) temp=temp*alph0 eps_error=max(abs(temp-u(i,j),eps_error) u(i,j)=temp enddo iteration=iteration+1enddoprint*,iteration,eps_errordo j=1,line,1 do i=1,mpoint,1 field(i,j)=u(i,j) enddoenddodeallocate(u,stat=ierr)endsubroutine!*2d网格数据无约束最小曲率迭代结束*!*2d余弦扩边开始*subrout

41、ine expan_cos2(m1,m2,m3,m4,n1,n2,n3,n4,g,flag) parameter (pi=3.141592654) integer flag integer m1,m2,m3,m4 integer n1,n2,n3,n4 real g(m1:m4,n1:n4) real g1(m1:m4,n1:n4) real g2(m1:m4,n1:n4) real sum integer num sum=0.0 num=0 if(flag=1) then !扩边端点值取边界平均值 do i=m2,m3 sum=sum+g(i,n2)+g(i,n3) num=num+2 en

42、d do do i=n2+1,n3-1 sum=sum+g(m2,i)+g(m3,i) num=num+2 end do do i=m1,m4 g(i,n1)=sum/num g(i,n4)=sum/num end do do i=n1+1,n4-1 g(m1,i)=sum/num g(m4,i)=sum/num end do else if(flag=2) then !扩边端点值取全部数据平均值 do i=m2,m3 do j=n2,n3 sum=sum+g(i,j) num=num+1 end do end do do i=m1,m4 g(i,n1)=sum/num g(i,n4)=sum

43、/num end do do i=n1+1,n4-1 g(m1,i)=sum/num g(m4,i)=sum/num end do else if(flag=3) then !扩边端点值取常数， do i=m1,m4 g(i,n1)=0 g(i,n4)=0 end do do i=n1+1,n4-1 g(m1,i)=0 g(m4,i)=0 end do end if!扩边 g1=g g2=g do j=n2,n3 do i=m1+1,m2-1 g1(i,j)=(g1(m2,j)-g1(m1,j)*cos(m2-i)*1.0/(m2-m1)*pi/2.0)+g1(m1,j) end do do

44、i=m3+1,m4-1 g1(i,j)=(g1(m4,j)-g1(m3,j)*cos(m4-i)*1.0/(m4-m3)*pi/2.0)+g1(m3,j) end do end do do i=m1,m4 do j=n1+1,n2-1 g1(i,j)=(g1(i,n2)-g1(i,n1)*cos(n2-j)*1.0/(n2-n1)*pi/2.0)+g1(i,n1) end do do j=n3+1,n4-1 g1(i,j)=(g1(i,n4)-g1(i,n3)*cos(n4-j)*1.0/(n4-n3)*pi/2.0)+g1(i,n3) end do end do do i=m2,m3 do

45、j=n1+1,n2-1 g2(i,j)=(g2(i,n2)-g2(i,n1)*cos(n2-j)*1.0/(n2-n1)*pi/2.0)+g2(i,n1) end do do j=n3+1,n4-1 g2(i,j)=(g2(i,n4)-g2(i,n3)*cos(n4-j)*1.0/(n4-n3)*pi/2.0)+g2(i,n3) end do end do do j=n1,n4 do i=m1+1,m2-1 g2(i,j)=(g2(m2,j)-g2(m1,j)*cos(m2-i)*1.0/(m2-m1)*pi/2.0)+g2(m1,j) end do do i=m3+1,m4-1 g2(i,j

46、)=(g2(m4,j)-g2(m3,j)*cos(m4-i)*1.0/(m4-m3)*pi/2.0)+g2(m3,j) end do end do g=(g1+g2)/2.0end subroutine expan_cos2!*2d余弦扩边开始*!*扩边子程序集结束*!*计算导数因子开始*subroutine deriv_operator_sub(dx,dy,m,n,deriv_operator_x,deriv_operator_y,deriv_operator_z) implicit nonereal dx,dyreal dzinteger m,n,i,jreal w(1:m,1:n),u(1:m),v(1:n)real deriv_operator_z(1:m,1:n)real deriv_operator_x(1:m,1:n)real deriv_operator_y(1:m,1:n)call wave2d_sub(u,v,w,m,n,dx,dy) !计算垂向导数因子,实数deriv_operator_z=w !计算x和y导数因子，虚数 do j=1,n do i=1,m deriv_operator_x(i,j)=u(i)