计算最优化离散波数的优化算法

第27卷第1期物探化探计算技术2005年2月文章编号: 1001- 1749(2005)01- -0034- 05计算最优化离散波数的优化算法柳建新，刘海飞(中南大学信息物理工程学院，长沙410083)摘要:在阐述用最优化方法计算离散波数的基础上,对波数初值的给定及偏导数矩阵的计算方法作了进-步的改进,使其在数学推理上更加严密;在计算量和计算精度方面也有了很大的改善,并简化了程序设计。通过试算发现,反付氏变换的均方误差随波数个数的变化规律,即在变化曲线上存在转折点,在转折点之后误差趋于平稳变化。选取转折点处的波数个数作为反付氏变换的波数个数,这样在正演、反演过程中,既保证了计算精度,又节约了计算时间。关键词:离散波数;偏导数矩阵;均方误差;正演;反演中图分类号: O241. 1文献标识码: A是文献[3]没有说明波数初值的给定方法。付氏电前言位对波数的偏导数也是用差商形式计算的,那么差商的步长选取就有了一定的人为性。对此,作者对用有限元法解点源二维问题时,要用到反付氏上述不足加以改进,如用计算等比离散波数的方法变换,用数值方法进行反付氏变换时,波数λ的选给定波数的初值,用解析方法计算付氏电位对波数取是保证计算精度和节约计算时间的主要问题。罗的偏导数,使其在数学推理上变的更加严密;并且延钟[.2在进行反付氏变换时,根据零阶修正贝塞改善计算量和计算精度,简化程序设计。通过试算，尔函数K。(x)曲线的特点,构造线性函数和负指数发现了反付氏变换的均方误差随波数个数的变化函数,然后对其进行分段积分来完成反付氏变换。的规律,即在变化曲线上存在转折点，在转折点之该方法在反付氏变换过程中选择的波数呈等比数后误差趋于平稳变化。选取转折点处的波数个数作列,通常称此波数序列为等比离散波数。它的缺点为反付氏变换的波数个数,这样在正演、反演过程就在于计算波数时仅考虑了最大和最小电极距,只中,既保证了计算精度,又节约了计算时间。要最大和最小电极间距不变,计算的波数序列也就不变。所以用等比离散波数进行反付氏变换时,需计算最优化离散波数的基本原理要较多的波数个数才能满足模拟精度,那么正演模拟的计算量也将随着波数个数的增加而增加。为根据文献[3]给出的用最优化方法计算离散波此，徐世浙[8]利用最优化方法计算离散波数,它的数的基本原理,U(x,y,z)表示点源二维介质电场优点就在于将电极距序列的每个电极距都参与计中的三维电位,V(x,),z)表示沿走向方向进行付算,最终使反付氏变换的计算结果达到最优,用该中国煤化工正演是通过先求出二维方法计算出来的离散波数被称为最优化离散波数。YHcNMHGt氏变换.最终求得实际该方法可以在较少的波数个数下就能达到较高的的三维电位U(x,y,z)来实现的。在主剖面上,当y模拟精度,它是一种很有效的反付氏变换方法。但=0 时,反付氏变换为基金项目:“开方姻屠科技攻关项目资助(2001BA609A-06)收稿日期: 2004- 06- 031期柳建新等:计算最优化离散波数的优化算法U(x,0,z)=三V(x,n,z)di(1)的最优化离散波数入和反付氏变换系数g;代入式.πJo(2),进行反付氏变换。式(1)的积分可以写成U(r)= 2 V(r,x)g;+g.+1(2)2最优化离散波数的计算过程其中x(j=1,2..,n)是 最优化离散波数; g;(j计算最优化离散波数可分以下三个步骤进行=1,2,..n+1)是反付氏变换系数; r是主剖面上2.1给定波 数初值测点到供电点的距离。用求等比离散波数的方法给定一组初始波数研究的目的是通过有效的方法选择合适的入;λ。根据文献[1]选择反付氏变换的积分界限[a,和g,使式(2)在，的一定范围内尽可能准确。但是6],通常积分下限a=0. 1,积分上限b=3.0。设在在一般情况下,函数U、V的表达式是未知的,因此已给定模拟尺度条件下,所用电极装置的最大和最无法直接利用式(2)求),和gj。由此利用均匀半空小电极距分别为rma和rmin,则最小波数λ和最大.间模型的反付氏变换公式代替式(1),有波数A可分别选为bK。(2r)dλ(3)Tmin其中K。是第二类零阶修正贝塞尔函数。将式(3)接着计算中间段的离散波数,为消除文献[1,写成形如式(2)的形式,有2]在计算等比离散波数时,由于先给定常数c后计算波数个数n,将对波数个数带来较大的舍入误Ko(rA;)g;+gn+1(4)差[4],故文献[4]先给定波数个数n,然后再根据式为了在不同的r下,有相近的相对误差,将式(10)计算常数cIn),- Inλ,(4)写成(10)n- 11= jrKo(rA;)g;+rgn+1(5)把常数c代入式(11),进而可以得到中间的等比离然后给定一-组电极距序列r;(i=1,2...,m),将(5)散波数式右端写成矩阵和向量形式,并令其等于向量V,λ;+1=λe° j=1,2,.,n-2得2.2 计算反付氏变换系数AG=V(6)计算出波数后,接着再求出矩阵A。由目标函其中A= (a;)mx(n+1),为一mX(n+1)阶矩阵。当数φ对g;取极小,并令其等于零,即j=1,2...n和i=1,2，...m时，元素aj=r,Ko(r,A);当j=n+1和i= 12...m时,元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计算最优化离散波数式(7)中的符号含义同上。由于目标函数φ受λ,和g;两个因素的共同决构造目标函数φ中国煤化工-组波数入;下,目标函CNMHGv在一组初始2(处展.φ=(I-V)"(I-V)=成泰勒级数,并取入的一次项,有(I- AG)"(I- AG)(8)V=V。+aV .(13)分别选取方科据使目标函数中达到极小,其中1为单位列向量。当目标函数φ达到极小时,把得到将式(13)代入到目标函数φ中,有●36●物探化探计算技术27卷φ=(I-V。V8xr根据贝塞尔函数的递推公式[5]"K,(x)]=-x "K+1(x)(18)(I-V。-Vax)(14)此时的φ是δλ的函数,由φ的极小,可决定当n=0时,有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)个计算过程,直到迭代后的均方误差ε= V中/m小于事先给定的允许误差限为止,最终得到的波数序对于上述第二类修正贝塞尔函数K。(x)和K:列作为最优化离散波数。作者对上述超定方程组式(x),可采用下面两个近似公式进行计算,计算精(12)和式(15)采用奇异值分解算法进行求解,计算度较高0]。当x≥2时,有精度较高。Ko(x)=(1.25331414+(2/x)3偏导数矩阵Bmn的计算(一0.07832358 + ((2/x)(0. 02189568 +(2/x)(-0.01062446+ (2/x)(0.00587872+文献[3]在计算偏导数矩阵Bmxn 时,通过给定(2/x)(- 0.0025154+(2/x)-组波数.....,经22过程可计算出V。0.0005320)))))/(x●e*)再给定一组波数的扰动量,文献[3]取扰动量0x;=0.1入。再计算出另一组波数λ;+0),再经2.2过.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")不难看出，经两次2.2过程,才计算出付氏电当x<2时,有位对波数的差商形式,这无疑增加了计算量。由于扰动量是人为给出的,那么扰动量的大小决定着差K,(x)=-ln(x/2)(1.0+ (x/3.75)2商对偏导数的近似程度.扰动量不同,近似程度也(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期柳建新等:计算最优化离散波数的优化算法●37●(x/3.75)2(0. 00301532 + (x/3.75)25]区间内,误差曲线的变化较快,而在n∈[5,10]0. 003211)))))n(x/2) +(1/x)区间内,误差曲线变化甚微,基本是同一个数量级(1. 00+十(x/2)*(0.15443144+(x/2)2的计算精度。这说明在n=5时,是误差随波数个数(一0.67278579+ (.x/2)2(-0.18156897+变化的一个转折点,则n=5可作为最佳波数个数(x/2)2(-0. 01919402 + (x/2)2的选取值。(- 0. 00110404- (x/2)*0.00004686)))))5数值模拟4最佳波数个数的选择为验证本文算法的有效性,仍然以上述电极距序列为例,对半空间均匀介质模型用最优化离散波在有限元正演模拟过程中,波数个数的多少直数进行反付氏变换计算，并与文献[1,2]用等比离接影响着数值模拟的计算量和计算精度。为保证数散波数进行反付氏变换的结果进行对比,得到均匀值模拟的精度,又不要徒劳的增加计算量,所以,选介质模型数值模拟结果如下页表1。择一个合理波数个数是至关重要的。从表1中可以看出,在均匀介质条件下,用最在均匀半空间地电模型条件下,给定-组电极优化离散波数进行反付氏变换计算,较少的波数就距AO=1. 5、2. 5、4.6、9、15.25、40、65、100、150、可以达到较高的模拟精度。而用等比离散波数进行220、260、300,用公式U;=Iρ/2πr;计算每个电极距模拟计算,即使在选取较多的波数个数的情况下，的理论电位,不失一般性，设Ip/2π=1。采用最优计算精度依然不及选取较少的最优化离散波数个化离散波数计算每个电极距的近似电位U,并计数的计算结果,这个对比结果是非常令人满意的。算近似解的均方误差6结论(U,-U,)*/m(22)通过对计算最优化离散波数的几个关键性技i=1,2... ,m术问题的研究,并经过编程试算,可以得到以下几式中m 为电极距个数。如图1所示，观察均方误点结论。差随波数个数的变化曲线。在波数个数为n∈[2,(1)用计算等比离散波数的方法给定计算最优.5 [化离散波数的初始波数,可以避免因手工输入的初始波数偏离最优解过大，而导致计算过程收敛慢的影响。2.0 t(2)用解析方法计算付氏电位对波数的偏导数矩阵Bmxn ,使其在数学推理上更加严密,在计算量i1.5t和计算精度方面也有了很大的改善,同时也简化了程序设计。.0 t(3)用最优化离散波数进行反付氏变换时,选择波数个数N=5作为最佳波数个数,即使再增加0.5波数个数,模拟精度也无明显改善,但是计算量却明显增加。0.0 t中国煤化工比,用最优化离散波数YHCNMHG的波数个数就能达到较n/波数个数高的模拟精度,而且节省了正演模拟的计算量,这图1均方误差随波数个数的变化曲线将为长数据断面的二维正反演提供了可能。Fig. 1 The curve changed of square errors VS_ the, number of wave万万数据●38●物探化探计算技术27卷表1对均匀介质模型数值模拟结果Tab.1 Numerical results for the model of homogenized medium等比离散波数最优化离散波数极距理论值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] 南京工学院数学教研组编.数学物理方程与特殊函参考文献:数[M].北京:高等教育出版社,1994.[1]罗延钟.电子计算机在电法勘探中的应用[M].武汉[6]徐士良. C常用算法程序集(第二版)[M].北京:清地质学院出版社,1987..华大学出版社，1996.[2]罗延钟.关于用有限元法对二维构造作电阻率法模拟的几个问题[J].地球物理学报,1986 ,29(6):11.作者简介:柳建新(1962~),男，湖南岳阳人,教[3]徐世浙. 点电源二维电场问题中付氏变换的波数K授,博士生导师,中南大学信息物理工程学院副院的选择[J].物探化探计算技术,1988,10(3):8. .长，中国地球物理学会海洋专业委员会理事、湖南[4]刘海飞.高密度电阻率法数据处理方法研究[D].中省地球物理学会副理事长，长期从事电法、电磁法南大学硕士论文,2004.方面的理论与应用研究。中国煤化工MHCNMHG