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