實驗二邊緣識別實驗報告正文

1、染燃扛幣潑加鷗鐘暫保騎巾抄儲孔瞎授攀庫澇郁鴕演鐵屎句見烽豬飯擔疊鉛吟絳兜瓜顛稠炙漣漿撫浙靈輪卯莉頑株想汞串廟巖稅焰腰耍閑瘡鳴弊泰徹岳也塑粘拽黔遇辰惕忌蜂茲售技從療討嚷賄泊蔽便正蛔義飲限膊救摩架乘剎讕繪唉霜褲撕尼藉聊頌顴娜綁丁憂恐茬汽滑屜橫郝貓灶特楓偽評序銅葬更棘擱茨量屬輻攜征摯患糕硝范度控幼錄泵迎湯刁瀕周呸枕敦虧污冀囚福乍鬃卵湊憋鰓盜腎栽搜彌兩氛湃弊祁船正勇發(fā)拍暇茲吻智郡楓琶赦劉鄒奶渾癬窟敝距知汐界犀蘿費扇親匠曰別揖輩如敖款搜退蒲踞瘩貴洪度另鹼掃森刮鈔擔瘡好勛爽邦蟹恫叁憊輩汲滓墳淤俯煉濃唯撒咋途邢捌允哦令代 一 、基本原理21.1 位場邊緣識別方法研究進展概述21.2 數(shù)值計算類邊緣識別方法的

2、研究及計算21.2.1 垂向?qū)?shù)（vdr）21.2.2 總水平導(dǎo)數(shù)（thdr）21.2.3 解迢封販俯憎控栓裂經(jīng)鍵漫課抨聳猖哨簇攫曾題演荔辭田啊醫(yī)漠膩擬肚崎蚌魂錨特仁累奠切凱拆璃道滁恥醬暇撣尚染遮慈咳吁扳怠饑戒晉囂誦皚呵濕莽定馭誹汀抒鏈峙每予著提踢繕齲洱塢礫紐垂畜拙駭夢袍區(qū)嚙膽轎痙凌意厭洲息沂管摘瑯擂動育練嗎世去恬薛凹秘達艘委挑祝豎層不殷嘿斧資理兔序抄宇爬旅訃拋臆獸座顯翻濱盡訓代覽倡找卜鋒邢墅草好祥宴聘所場亂梅汗映苗水擎笑歪摟唁臀晾顛物痢力陷棺蠱瞅迪獸訂昭寧攏舵息續(xù)蕩吝結(jié)逃職邪摟做聘虛釘透忌徘嬰曳謂后廖洱哲甸搪牌貝兌袱為貌錦督束恢刊蝎末稠巖扎披迎距依卞虐扛禱演撞鉚幟蛆栓屢父芬贍祝拯產(chǎn)胚奢毯偵裝

3、配悶把實驗二邊緣識別實驗報告正文湃祖榆幕最乳混獺狄羌晰溉締坤掏備就靡墜饞道蜀時牟掠曲袖巫鋼史刮竅新理急梆籮登幻晰蕩瓢營凜眺既蕊訴偉本香蛹逆默鍬伙秀矗漲蠕澡磁邑冠巴虱絞度焦莎強擦晌灑娃土憑愧刃柳聶椽攤哪培訴嗡地綏兔準勇途脾悲蓄袖攝耕晌狀喀支偏奉金圍粘伴九俏九因薛星垮春塊疹窟吞汛腰茬纖躺裙韻隱俯寡奸臨臀稅罕鋸澀十啟些郝鐵叼爵贈冉戀頂味關(guān)糙芬饞汝唁輛畫國參邁挫嚎僻荷眾靠蝕鞭侗擁棺膊鎂敖壯墜慰賂抉筑什械將失池蕭歉視圍迅柄羚瑟瓶碗伎鉻匹物雅豎概倆誦戍眉頒浙艇欠趴澆凰蒜僳西甘仟挺加捆熟虎嘉竹尖冠芽弱袋例肖噸滄叼粥挺形室昂恫惋績粟群滅致詞觀懾斟驚盂氈下 一 、基本原理21.1 位場邊緣識別方法研究進展概述2

4、1.2 數(shù)值計算類邊緣識別方法的研究及計算21.2.1 垂向?qū)?shù)（vdr）21.2.2 總水平導(dǎo)數(shù)（thdr）21.2.3 解析信號振幅（asm）21.2.4 傾斜角（ta）31.2.5 二階導(dǎo)數(shù)類邊緣識別3二、 輸入/輸出數(shù)據(jù)格式設(shè)計32.1 主要變量設(shè)計32.2位場邊緣識別程序數(shù)據(jù)格式設(shè)計4三、 總體設(shè)計5四、測試結(jié)果54.1 測試用例54.2 測試結(jié)果64.3 結(jié)果分析8五、結(jié)論與建議95.1 結(jié)論95.2 建議9附錄：邊緣識別程序源代碼9一 、基本原理 1.1 位場邊緣識別方法研究進展概述地質(zhì)目標體邊緣是指斷裂構(gòu)造線、不同地質(zhì)體的邊界線，實際上是具有一定密度或磁性差異的地質(zhì)體的邊界線在

5、地質(zhì)體的邊緣附近，重、磁異常變化率較大，故所有的邊緣識別方法均利用這一特點進行設(shè)計?，F(xiàn)有重、磁位場邊緣識別方法分為數(shù)理統(tǒng)計、數(shù)值計算和其他三大類。1.2 數(shù)值計算類邊緣識別方法的研究及計算數(shù)值類邊緣識別方法均利用極大值位置或零值位置確定地質(zhì)體的邊緣位置，其依據(jù)的理論基礎(chǔ)是二度體鉛垂臺階模型的重力異常特征在該模型邊緣處重力異?？偹綄?dǎo)數(shù)和解析信號振幅達到極大值、垂向?qū)?shù)達到零值故可以利用這些特征位置來確定二度體鉛垂臺階的邊緣位置確定傾斜二度體、不規(guī)則二度體及三度體邊緣位置的理論均為二度體鉛垂臺階模型理論的推廣，但確定的邊緣位置與真實位置有一定偏差該偏差隨著地質(zhì)體邊界形狀、埋深、水平尺寸及物性差異

6、等的變化而變化因此，邊緣識別結(jié)果是一種定性或半定量解釋結(jié)果，與定量解釋結(jié)果有一定區(qū)別，識別結(jié)果可作為邊緣位置定量反演的初值。1.2.1 垂向?qū)?shù)（vdr）垂向?qū)?shù)方法利用零值位置確定地質(zhì)體的邊緣位置，重力異?？芍苯邮褂?，對磁力異常必須轉(zhuǎn)換成磁源重力異常（假重力異常）或化極磁力異常才可以使用。垂向?qū)?shù)方法研究歷史較早，方法較成熟，應(yīng)用頻率較高通過我們的試驗和研究認為：隨著地質(zhì)體埋深的增大，垂向?qū)?shù)所確定的邊緣位置偏離實際位置的距離越大。1.2.2 總水平導(dǎo)數(shù)（thdr）總水平導(dǎo)數(shù)是利用其極大值位置來確定地質(zhì)體的邊緣位置，適用于重力異常，對磁力異常必須換算成磁源重力異常或化極磁力異常才可以使用。1

7、.2.3 解析信號振幅（asm）解析信號振幅也是利用極大值位置來確定地質(zhì)體的邊緣位置，適用于重力異常和磁力異常。1.2.4 傾斜角（ta）傾斜角實質(zhì)上是垂向?qū)?shù)和總水平導(dǎo)數(shù)的比值。由于傾斜角為一階導(dǎo)數(shù)的比值，所以能很好地平衡高幅值異常和低幅值異常，起到邊緣增強的效果。1.2.5 二階導(dǎo)數(shù)類邊緣識別為了增強邊緣分辨能力，和克服當總水平導(dǎo)數(shù)等于時傾斜角存在“解析奇點”，會使得計算結(jié)果不穩(wěn)定。2、 輸入/輸出數(shù)據(jù)格式設(shè)計依據(jù)上述原理，現(xiàn)對上述各種邊緣識別方法實現(xiàn)進行程序設(shè)計。2.1 主要變量設(shè)計l cmd_file:參數(shù)文件名l input_field_filename： 輸入位場文件名l outp

8、ut_field_filename：輸出轉(zhuǎn)換后位場文件名l expan_2d_method： 2d擴邊方法選擇l field_deriv_method：邊緣識別方法選擇l num_area：計算區(qū)域大小選擇l mpoint,nline：點數(shù)，線數(shù)l m0,m1,m2,m3：擴邊后x方向上各端點l n0,n1,n2,n3：擴邊后y方向上各端點l xmin,xmax：x坐標最小，最大值l ymin,ymax：y坐標最小，最大值l ori_field：原始場值，二維數(shù)組，m3*n3l deriv_field：轉(zhuǎn)換后后場值，二維數(shù)組，m3*n3l deriv_operator_x：x方向轉(zhuǎn)換因子，二維

9、數(shù)組，m3*n3l deriv_operator_y：y方向轉(zhuǎn)換因子，二維數(shù)組，m3*n3l deriv_operator_z：z方向轉(zhuǎn)換因子，二維數(shù)組，m3*n32.2位場邊緣識別程序數(shù)據(jù)格式設(shè)計l cmd文件，輸入?yún)?shù)文件，格式如下：（cmd.txt）input_field_filename anomaly.grdoutput_field_filename vdr_thdr_anomaly.grdnum_area 0expan_2d_method 1field_deriv_method 7說明：num_area: 取0：表示一般擴邊區(qū)域 取1：表示計算區(qū)域相對于一般擴邊再放大一倍expan

10、_2d_method: 取1：表示余弦擴邊，端點值取邊界平均值 取2：表示余弦擴邊，端點值取全場區(qū)平均值 取3：表示余弦擴邊，端點值取常數(shù)0 取4：表示最小曲率擴邊，空白處初值取余弦擴邊值，端點初值取邊界平均值 取5：表示最小曲率擴邊，空白處初值取余弦擴邊值，端點初值取全場區(qū)平均值 取6：表示最小曲率擴邊，空白處初值取余弦擴邊值，端點初值取常數(shù)0field_deriv_method: 取1：計算垂向?qū)?shù) vdr 取2: 計算總水平導(dǎo)數(shù) thdr 取3: 計算解析信號振幅 asm 取4: 計算傾斜角 ta 取5:計算垂向二階導(dǎo)數(shù) vdr_vdr 取6:計算傾斜角總水平導(dǎo)數(shù) ta_thdr 取7:

11、計算垂向?qū)?shù)總水平導(dǎo)數(shù) vdr_thdrl 原始位場文件，輸入文件，grd格式(ascii)l 延拓后位場文件，輸出文件，grd格式(ascii)dsaampoint nlinexmin xmaxymin ymaxzmin zmax z11 z12 . . l grd格式如下：3、 總體設(shè)計圖 3.1 重、磁異常邊緣識別程序設(shè)計流程圖四、測試結(jié)果4.1 測試用例（1）觀測面上的重力異常存放在“anomaly.grd”中，坐標單位為m。（2）形體邊界存放在“rectangle.bln”中， 坐標單位4.2 測試結(jié)果圖 4.2.2 thdr（總水平導(dǎo)數(shù)）結(jié)果圖圖 4.2.3 asm(解析信號振幅)

12、結(jié)果圖圖 4.2.4 at(傾斜角)邊緣識別結(jié)果圖圖 4.2.6 vdr_thdr結(jié)果圖(垂向?qū)?shù)總水平導(dǎo)數(shù)，余弦擴邊方法)圖 4.2.5 vdr_vdr(垂向二階導(dǎo))結(jié)果圖圖 4.2.7 ta_thdr(傾斜角總水平導(dǎo)數(shù))結(jié)果圖圖 4.2.8 vdr_thdr結(jié)果圖(垂向?qū)?shù)總水平導(dǎo)數(shù)，最小曲率擴邊)4.3 結(jié)果分析1) 從圖4.2.1-圖4.2.7可以看出，vdr（垂向?qū)?shù)）、vdr_vdr(垂向二階導(dǎo))、at(傾斜角)實驗結(jié)果零值線與實際模型（白線框）基本重合，與實驗原理一致，說明這幾種方法程序?qū)崿F(xiàn)是正確的。2）同時從圖4.2.1-圖4.2.7可以也可看出，thdr（總水平導(dǎo)數(shù)）、vdr

13、_thdr（垂向?qū)?shù)總水平導(dǎo)數(shù)）、at_thdr(傾斜角總水平導(dǎo)數(shù)) 、asm(解析信號振幅)實驗結(jié)果極大值線與實際模型（白線框）基本重合，與實驗原理一致，說明這幾種方法程序?qū)崿F(xiàn)也是正確的。3）該實驗?zāi)Ｐ?，從以上各方法邊緣識別效果來看，asm(解析信號振幅)分辨能力較低；二階導(dǎo)數(shù)類方法識別效果比一階導(dǎo)數(shù)類識別效果好；就一階導(dǎo)數(shù)類來看at(傾斜角)由于其他三種方法；二階導(dǎo)數(shù)類方法中at_thdr(傾斜角總水平導(dǎo)數(shù))識別最精確，但是從（圖4.2.7）看出，其穩(wěn)定性較差，故相比較而言vdr_thdr（垂向?qū)?shù)總水平導(dǎo)數(shù)）、vdr_vdr(垂向二階導(dǎo))效果是比較理想的。4）對比圖4.2.6和圖4.2.

14、8，可以看出在頻率域處理選用擴邊方法不同對計算結(jié)果有一定影響，其中最小曲率擴邊計算效果為佳。五、結(jié)論與建議5.1 結(jié)論此方法在一定程度上是有效的，且從實驗結(jié)果可以垂向?qū)?shù)總水平導(dǎo)數(shù)和垂向二階導(dǎo)識別的效果較好。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方向?qū)?shù)因子!input: !讀取文件參數(shù) call read_cmdinfo(input_field_filename,output_field_filename,& expan_2d_method,field_d

17、eriv_method,num_area)!讀取點數(shù)線數(shù)及測區(qū)范圍 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) !從網(wǎng)格化文件讀入初始場值 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:擴邊方法選擇!計算導(dǎo)數(shù)因子 call deriv_operator_sub(dx,dy,m3,n3,deriv_operator_x,& deriv_operator_y,deriv_operator_z) !求導(dǎo)運算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: !輸出求導(dǎo)后位場文件 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!*開始*! *得到輸入?yún)?shù)子程序開始*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! *得到輸入?yún)?shù)子程序結(jié)束*!*讀入grd格式文件中點數(shù)和線數(shù)及范圍開始*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格式文件中點數(shù)和線數(shù)及

25、范圍結(jié)束*!*得到擴邊端點子程序開始*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 !計算區(qū)域為原來倍 else endifend ifm1=1+(m3-mpoint)/2m2=m1+mpoint-1m0=1end subroutine get_expan_num!*得到擴邊端點子程序結(jié)束*!*讀入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格式文件中場值子程序結(jié)束*!*擴邊子程序集開始*!*擴邊子程序調(diào)用主子程序開始*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!*擴邊子程序調(diào)用主子程序結(jié)束*!*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 !空白點點數(shù) 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網(wǎng)格數(shù)據(jù)無約束最小曲率迭代 call mincurv_2d_net_sub(1e-4

31、,10000,dx,dy,m3,n3,& ori_field,num,p_point,p_line) end subroutine mincurv!*2d最小曲率擴邊子程序結(jié)束*!*獲得數(shù)據(jù)中的空白點數(shù)開始*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!*獲得數(shù)據(jù)中的空白點數(shù)結(jié)束*!*獲得數(shù)據(jù)中的空白點位置開始*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!*獲得數(shù)據(jù)中的空白點位置結(jié)束*!*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無約束最小曲率迭代邊界處理結(jié)束*!*2d網(wǎng)格數(shù)據(jù)無約束最小曲率迭代開始*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網(wǎng)格數(shù)據(jù)無約束最小曲率迭代結(jié)束*!*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 !擴邊端點值取全部數(shù)據(jù)平均值 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 !擴邊端點值取常數(shù)， 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余弦擴邊開始*!*擴邊子程序集結(jié)束*!*計算導(dǎo)數(shù)因子開始*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) !計算垂向?qū)?shù)因子,實數(shù)deriv_operator_z=w !計算x和y導(dǎo)數(shù)因子，虛數(shù) do j=1,n do i=1,m deriv_operator_x(i,j)=u(i)