計算最優(yōu)化離散波數的優(yōu)化算法

第27卷第1期物探化探計算技術(shù)2005年2月文章編號: 1001- 1749(2005)01- -0034- 05計算最優(yōu)化離散波數的優(yōu)化算法柳建新，劉海飛(中南大學(xué)信息物理工程學(xué)院，長(cháng)沙410083)摘要:在闡述用最優(yōu)化方法計算離散波數的基礎上,對波數初值的給定及偏導數矩陣的計算方法作了進(jìn)-步的改進(jìn),使其在數學(xué)推理上更加嚴密;在計算量和計算精度方面也有了很大的改善,并簡(jiǎn)化了程序設計。通過(guò)試算發(fā)現,反付氏變換的均方誤差隨波數個(gè)數的變化規律,即在變化曲線(xiàn)上存在轉折點(diǎn),在轉折點(diǎn)之后誤差趨于平穩變化。選取轉折點(diǎn)處的波數個(gè)數作為反付氏變換的波數個(gè)數,這樣在正演、反演過(guò)程中,既保證了計算精度,又節約了計算時(shí)間。關(guān)鍵詞:離散波數;偏導數矩陣;均方誤差;正演;反演中圖分類(lèi)號: O241. 1文獻標識碼: A是文獻[3]沒(méi)有說(shuō)明波數初值的給定方法。付氏電前言位對波數的偏導數也是用差商形式計算的,那么差商的步長(cháng)選取就有了一定的人為性。對此,作者對用有限元法解點(diǎn)源二維問(wèn)題時(shí),要用到反付氏上述不足加以改進(jìn),如用計算等比離散波數的方法變換,用數值方法進(jìn)行反付氏變換時(shí),波數λ的選給定波數的初值,用解析方法計算付氏電位對波數取是保證計算精度和節約計算時(shí)間的主要問(wèn)題。羅的偏導數,使其在數學(xué)推理上變的更加嚴密;并且延鐘[.2在進(jìn)行反付氏變換時(shí),根據零階修正貝塞改善計算量和計算精度,簡(jiǎn)化程序設計。通過(guò)試算，爾函數K。(x)曲線(xiàn)的特點(diǎn),構造線(xiàn)性函數和負指數發(fā)現了反付氏變換的均方誤差隨波數個(gè)數的變化函數,然后對其進(jìn)行分段積分來(lái)完成反付氏變換。的規律,即在變化曲線(xiàn)上存在轉折點(diǎn)，在轉折點(diǎn)之該方法在反付氏變換過(guò)程中選擇的波數呈等比數后誤差趨于平穩變化。選取轉折點(diǎn)處的波數個(gè)數作列,通常稱(chēng)此波數序列為等比離散波數。它的缺點(diǎn)為反付氏變換的波數個(gè)數,這樣在正演、反演過(guò)程就在于計算波數時(shí)僅考慮了最大和最小電極距,只中,既保證了計算精度,又節約了計算時(shí)間。要最大和最小電極間距不變,計算的波數序列也就不變。所以用等比離散波數進(jìn)行反付氏變換時(shí),需計算最優(yōu)化離散波數的基本原理要較多的波數個(gè)數才能滿(mǎn)足模擬精度,那么正演模擬的計算量也將隨著(zhù)波數個(gè)數的增加而增加。為根據文獻[3]給出的用最優(yōu)化方法計算離散波此，徐世浙[8]利用最優(yōu)化方法計算離散波數,它的數的基本原理,U(x,y,z)表示點(diǎn)源二維介質(zhì)電場(chǎng)優(yōu)點(diǎn)就在于將電極距序列的每個(gè)電極距都參與計中的三維電位,V(x,),z)表示沿走向方向進(jìn)行付算,最終使反付氏變換的計算結果達到最優(yōu),用該中國煤化工正演是通過(guò)先求出二維方法計算出來(lái)的離散波數被稱(chēng)為最優(yōu)化離散波數。YHcNMHGt氏變換.最終求得實(shí)際該方法可以在較少的波數個(gè)數下就能達到較高的的三維電位U(x,y,z)來(lái)實(shí)現的。在主剖面上,當y模擬精度,它是一種很有效的反付氏變換方法。但=0 時(shí),反付氏變換為基金項目:“開(kāi)方姻屠科技攻關(guān)項目資助(2001BA609A-06)收稿日期: 2004- 06- 031期柳建新等:計算最優(yōu)化離散波數的優(yōu)化算法U(x,0,z)=三V(x,n,z)di(1)的最優(yōu)化離散波數入和反付氏變換系數g;代入式.πJo(2),進(jìn)行反付氏變換。式(1)的積分可以寫(xiě)成U(r)= 2 V(r,x)g;+g.+1(2)2最優(yōu)化離散波數的計算過(guò)程其中x(j=1,2..,n)是 最優(yōu)化離散波數; g;(j計算最優(yōu)化離散波數可分以下三個(gè)步驟進(jìn)行=1,2,..n+1)是反付氏變換系數; r是主剖面上2.1給定波 數初值測點(diǎn)到供電點(diǎn)的距離。用求等比離散波數的方法給定一組初始波數研究的目的是通過(guò)有效的方法選擇合適的入;λ。根據文獻[1]選擇反付氏變換的積分界限[a,和g,使式(2)在，的一定范圍內盡可能準確。但是6],通常積分下限a=0. 1,積分上限b=3.0。設在在一般情況下,函數U、V的表達式是未知的,因此已給定模擬尺度條件下,所用電極裝置的最大和最無(wú)法直接利用式(2)求),和gj。由此利用均勻半空小電極距分別為rma和rmin,則最小波數λ和最大.間模型的反付氏變換公式代替式(1),有波數A可分別選為bK。(2r)dλ(3)Tmin其中K。是第二類(lèi)零階修正貝塞爾函數。將式(3)接著(zhù)計算中間段的離散波數,為消除文獻[1,寫(xiě)成形如式(2)的形式,有2]在計算等比離散波數時(shí),由于先給定常數c后計算波數個(gè)數n,將對波數個(gè)數帶來(lái)較大的舍入誤Ko(rA;)g;+gn+1(4)差[4],故文獻[4]先給定波數個(gè)數n,然后再根據式為了在不同的r下,有相近的相對誤差,將式(10)計算常數cIn),- Inλ,(4)寫(xiě)成(10)n- 11= jrKo(rA;)g;+rgn+1(5)把常數c代入式(11),進(jìn)而可以得到中間的等比離然后給定一-組電極距序列r;(i=1,2...,m),將(5)散波數式右端寫(xiě)成矩陣和向量形式,并令其等于向量V,λ;+1=λe° j=1,2,.,n-2得2.2 計算反付氏變換系數AG=V(6)計算出波數后,接著(zhù)再求出矩陣A。由目標函其中A= (a;)mx(n+1),為一mX(n+1)階矩陣。當數φ對g;取極小,并令其等于零,即j=1,2...n和i=1,2，...m時(shí)，元素aj=r,Ko(r,A);當j=n+1和i= 12...m時(shí),元dφ.--2A"(I- AG)=0素a;=r,o向量G= (gi+g...go.gn+1)"。向量V=有AG=I(12)V(V.,V,.,V,...v..)“,V;可表示為成立,符號含義同前。對式(12)求解,就可得到一組V;= Zangj(7)g;,根據式(7)可求出付氏電位Vi。=12.3計算最優(yōu)化離散波數式(7)中的符號含義同上。由于目標函數φ受λ,和g;兩個(gè)因素的共同決構造目標函數φ中國煤化工-組波數入;下,目標函CNMHGv在一組初始2(處展.φ=(I-V)"(I-V)=成泰勒級數,并取入的一次項,有(I- AG)"(I- AG)(8)V=V。+aV .(13)分別選取方科據使目標函數中達到極小,其中1為單位列向量。當目標函數φ達到極小時(shí),把得到將式(13)代入到目標函數φ中,有●36●物探化探計算技術(shù)27卷φ=(I-V。V8xr根據貝塞爾函數的遞推公式[5]"K,(x)]=-x "K+1(x)(18)(I-V。-Vax)(14)此時(shí)的φ是δλ的函數,由φ的極小,可決定當n=0時(shí),有8。為此,求φ對8λ的求導,并令其等于零,得到dK。(x)=-K;(x)(19)dxBmxn8hnx1=Cmx1(15)成立,所以有其中Bmxn= |JV; Cmx1=(I- Vo)mx1aK。(r;. ))=-r,K:(r;.入)(20))，從式(15)中解出δλ,于是得到一組新的29)成立,從而得到= λ{°十δλ .(16) .=-將.KK1(r;. A)●g;再以2()作為初始值,重復2.2章節和2.3章節兩(i=1,2,..m; j=1,2,.,n)個(gè)計算過(guò)程,直到迭代后的均方誤差ε= V中/m小于事先給定的允許誤差限為止,最終得到的波數序對于上述第二類(lèi)修正貝塞爾函數K。(x)和K:列作為最優(yōu)化離散波數。作者對上述超定方程組式(x),可采用下面兩個(gè)近似公式進(jìn)行計算,計算精(12)和式(15)采用奇異值分解算法進(jìn)行求解,計算度較高0]。當x≥2時(shí),有精度較高。Ko(x)=(1.25331414+(2/x)3偏導數矩陣Bmn的計算(一0.07832358 + ((2/x)(0. 02189568 +(2/x)(-0.01062446+ (2/x)(0.00587872+文獻[3]在計算偏導數矩陣Bmxn 時(shí),通過(guò)給定(2/x)(- 0.0025154+(2/x)-組波數.....,經(jīng)22過(guò)程可計算出V。0.0005320)))))/(x●e*)再給定一組波數的擾動(dòng)量,文獻[3]取擾動(dòng)量0x;=0.1入。再計算出另一組波數λ;+0),再經(jīng)2.2過(guò).K,(x)=(1.25331414+(2/x)程可計算出V,于是用下面的差商形式代替付氏(0. 23498619+ (2/x)(- 0. 03655620+電位對波數的偏導數.(2/x)(0. 01504268+ (2/x)(-0.00780353+(2/x)(0. 00325614- (2/x)0.00068245)))))/aV;_ V;-V;示Ox;(x●e")不難看出，經(jīng)兩次2.2過(guò)程,才計算出付氏電當x<2時(shí),有位對波數的差商形式,這無(wú)疑增加了計算量。由于擾動(dòng)量是人為給出的,那么擾動(dòng)量的大小決定著(zhù)差K,(x)=-ln(x/2)(1.0+ (x/3.75)2商對偏導數的近似程度.擾動(dòng)量不同,近似程度也(3.5156229+ (x/3.75)=(3. 0899424 +不同。對此，下面給出付氏電位對波數的偏導數的(x/3.75)2(1. 2067492+ (x/3.75)2解析算法。(0.2659732 + (x/3. 75)"(0. 0360768+由于V;是入;的函數,即V;=V;(...,x,.,(x/3. 75)*0. 0045813)))))- 0.57721566+A,), 在式(7)兩端分別對入;求導,得到中國煤化工2)*(0.23069756+MHCNMH G2)<(0.00262698+aV;_ a(a;8++..+a8;+..+ a;(n+1)8n+1)_(x72)2(0. 0001075+ (x/2)0.00007)))))V;a[r. Ko(r. );). gj]_K:(x)=x(0. 50+ (x/3.75)2 (0.87890594 +a;(x/3.75)2(0. 51498869+ (x/3.75)2，2) .●gj(17)(O.15084934 + (x/3. 75)2 (O.02658733 +1期柳建新等:計算最優(yōu)化離散波數的優(yōu)化算法●37●(x/3.75)2(0. 00301532 + (x/3.75)25]區間內,誤差曲線(xiàn)的變化較快,而在n∈[5,10]0. 003211)))))n(x/2) +(1/x)區間內,誤差曲線(xiàn)變化甚微,基本是同一個(gè)數量級(1. 00+十(x/2)*(0.15443144+(x/2)2的計算精度。這說(shuō)明在n=5時(shí),是誤差隨波數個(gè)數(一0.67278579+ (.x/2)2(-0.18156897+變化的一個(gè)轉折點(diǎn),則n=5可作為最佳波數個(gè)數(x/2)2(-0. 01919402 + (x/2)2的選取值。(- 0. 00110404- (x/2)*0.00004686)))))5數值模擬4最佳波數個(gè)數的選擇為驗證本文算法的有效性,仍然以上述電極距序列為例,對半空間均勻介質(zhì)模型用最優(yōu)化離散波在有限元正演模擬過(guò)程中,波數個(gè)數的多少直數進(jìn)行反付氏變換計算，并與文獻[1,2]用等比離接影響著(zhù)數值模擬的計算量和計算精度。為保證數散波數進(jìn)行反付氏變換的結果進(jìn)行對比,得到均勻值模擬的精度,又不要徒勞的增加計算量,所以,選介質(zhì)模型數值模擬結果如下頁(yè)表1。擇一個(gè)合理波數個(gè)數是至關(guān)重要的。從表1中可以看出,在均勻介質(zhì)條件下,用最在均勻半空間地電模型條件下,給定-組電極優(yōu)化離散波數進(jìn)行反付氏變換計算,較少的波數就距AO=1. 5、2. 5、4.6、9、15.25、40、65、100、150、可以達到較高的模擬精度。而用等比離散波數進(jìn)行220、260、300,用公式U;=Iρ/2πr;計算每個(gè)電極距模擬計算,即使在選取較多的波數個(gè)數的情況下，的理論電位,不失一般性，設Ip/2π=1。采用最優(yōu)計算精度依然不及選取較少的最優(yōu)化離散波數個(gè)化離散波數計算每個(gè)電極距的近似電位U,并計數的計算結果,這個(gè)對比結果是非常令人滿(mǎn)意的。算近似解的均方誤差6結論(U,-U,)*/m(22)通過(guò)對計算最優(yōu)化離散波數的幾個(gè)關(guān)鍵性技i=1,2... ,m術(shù)問(wèn)題的研究,并經(jīng)過(guò)編程試算,可以得到以下幾式中m 為電極距個(gè)數。如圖1所示，觀(guān)察均方誤點(diǎn)結論。差隨波數個(gè)數的變化曲線(xiàn)。在波數個(gè)數為n∈[2,(1)用計算等比離散波數的方法給定計算最優(yōu).5 [化離散波數的初始波數,可以避免因手工輸入的初始波數偏離最優(yōu)解過(guò)大，而導致計算過(guò)程收斂慢的影響。2.0 t(2)用解析方法計算付氏電位對波數的偏導數矩陣Bmxn ,使其在數學(xué)推理上更加嚴密,在計算量i1.5t和計算精度方面也有了很大的改善,同時(shí)也簡(jiǎn)化了程序設計。.0 t(3)用最優(yōu)化離散波數進(jìn)行反付氏變換時(shí),選擇波數個(gè)數N=5作為最佳波數個(gè)數,即使再增加0.5波數個(gè)數,模擬精度也無(wú)明顯改善,但是計算量卻明顯增加。0.0 t中國煤化工比,用最優(yōu)化離散波數YHCNMHG的波數個(gè)數就能達到較n/波數個(gè)數高的模擬精度,而且節省了正演模擬的計算量,這圖1均方誤差隨波數個(gè)數的變化曲線(xiàn)將為長(cháng)數據斷面的二維正反演提供了可能。Fig. 1 The curve changed of square errors VS_ the, number of wave萬(wàn)萬(wàn)數據●38●物探化探計算技術(shù)27卷表1對均勻介質(zhì)模型數值模擬結果Tab.1 Numerical results for the model of homogenized medium等比離散波數最優(yōu)化離散波數極距理論值N=9N=3N=4N=5N=61.50. 666670. 697190. 6546230. 6651730. 6669092.50.40. 415960. 4107860. 4024690. 399460. 39997140.250. 258020. 2507770. 2484820. 2504090. 25016860. 1666670. 171340. 1631290. 166441 0. 1667520. 16640490. 111110. 114050. 1099560. 111810. 1107780. 11120910. 0666660. 06840. 06821020. 066487 0. 06682630. 0667302250.04 .0. 041030. 04032630. 0398521| 0. 04003580. 0399359 .4(0. 0250. 025640. 02443470. 0251403 0. 02490690. 0250341650. 0153840. 015770. 01523470. 0153787 0. 01541990. 015370.0111110. 011390. 01124690. 011067| 0. 01113630. 0110991200. 0083330. 008540. 00847080. 0083275| 0. 00832170. 00834261500. 0066660. 006830. 00670540. 0066829| 0. 00665070. 00667751800. 0055550. 005690. 00551510. 0055688 0. 00555130. 00555742200. 0045450. 004650. 00447180. 0045426 0. 00455430. 00453792600. 0038460. 003930. 00381070. 0038369 0. 00385380. 00383993000. 0033330. 00340. 00338290. 0033366| 0. 00332740. 003386均方誤差(%)0. 89670.4173730. 08384470. 0205270. 008577[5] 南京工學(xué)院數學(xué)教研組編.數學(xué)物理方程與特殊函參考文獻:數[M].北京:高等教育出版社,1994.[1]羅延鐘.電子計算機在電法勘探中的應用[M].武漢[6]徐士良. C常用算法程序集(第二版)[M].北京:清地質(zhì)學(xué)院出版社,1987..華大學(xué)出版社，1996.[2]羅延鐘.關(guān)于用有限元法對二維構造作電阻率法模擬的幾個(gè)問(wèn)題[J].地球物理學(xué)報,1986 ,29(6):11.作者簡(jiǎn)介:柳建新(1962~),男，湖南岳陽(yáng)人,教[3]徐世浙. 點(diǎn)電源二維電場(chǎng)問(wèn)題中付氏變換的波數K授,博士生導師,中南大學(xué)信息物理工程學(xué)院副院的選擇[J].物探化探計算技術(shù),1988,10(3):8. .長(cháng)，中國地球物理學(xué)會(huì )海洋專(zhuān)業(yè)委員會(huì )理事、湖南[4]劉海飛.高密度電阻率法數據處理方法研究[D].中省地球物理學(xué)會(huì )副理事長(cháng)，長(cháng)期從事電法、電磁法南大學(xué)碩士論文,2004.方面的理論與應用研究。中國煤化工MHCNMHG