各向异性弹性参数的广义非线性反演方法.pdf
第30卷第3期2007年6月勘探地球物理进展Progress in Exploration GeophysicsVol.30,No.3Jun.,2007收稿日期:2006208216;改回日期:2006212231。第一作者简介:孙银行(1978),男,中国石油大学(华东)地球资源与信息学院地球探测与信息技术专业在读硕士,主要从事地球物理各向异性方面的研究。文章编号:167128585(2007)0320179206各向异性弹性参数的广义非线性反演方法孙银行,乐友喜(中国石油大学(华东)地球资源与信息学院,山东东营257061)摘要:在非线性最优化理论的基础上结合地球物理反演问题,提出各向异性弹性参数的广义非线性反演方法。结合各向异性弹性参数反演,推导出反演所需的梯度向量,梯度向量中Jacobi矩阵各个元素的求取是利用与正演完全相同的方法,并构造了一个公式,从单位阵开始来修改函数的海色矩阵的逆矩阵,这样,不用求海色矩阵及其逆阵就能求函数的极小值,减少了计算量,给出了实际的反演结果。地质模型的反演结果证明,该方法是可行的,且反演精度高,具有较好的发展前景。关键词:各向异性;弹性参数;全波场反演;海色矩阵中图分类号:P631.443文献标识码:A 地面叠前全波场炮集记录是地震勘探中采集的第一手原始资料,它包含面波、直达波、各类反射波、衍射波及其他一切波动效应,含有最充分的地震波运动学和动力学信息,因此,采用叠前全波场资料来反演各向异性弹性参数,对于减少反演的多解性、提高解的精度有着重要意义。并且,叠前全波场反演无需进行波场的识别与分离及其他有关处理工作,因而极大地减轻了处理中的困难和处理人员的工作量,也减少了处理中的人为因素。叠前全波场反演的优越性已引起了研究者的兴趣,杨顶辉等1用有限差分法对全波场各向异性弹性参数反演进行了研究,但只实现了用qSH波资料反演一个2层TI模型的垂向和水平速度;张秉铭等2在杨顶辉等的研究基础上做了进一步改进,利用两分量炮集记录反演了一个3层TI模型的C11和C33两个弹性参数。张美根等3用有限元法对全波场各向异性弹性参数反演进行了研究,反演了TI介质与qP波和qSV波有关的C11,C13,C33和C55四个弹性参数。以上工作在反演时用的都是广义线性反演。广义线性反演将非线性的地球物理问题用泰勒公式展开,去掉二次以上的项,得到线性化的系统,然后采用迭代法求解。该方法简单易行,计算效率非常高,因而长期以来深受人们的青睐。经过多年的研究,该方法在解的惟一性和稳定性等方面取得了实质性的进展,在许多地球物理反演问题中取得了较好的应用成果。随着地球物理学的发展,人们发现地球物理问题往往是非线性的,去掉二次以上的项所得到的线性化系统与实际问题有较大的偏差,在进行反演的过程中会出现一些难以解决的问题,例如反演时阻尼系数难以确定,收敛速度很慢,陷入局部极小的可能性增大,解估计的分辨率难以有充分的保证等等。如果在泰勒展式中取二次项得到一个二阶非线性系统,导出广义非线性反演方法,那么其解空间的性质、状态均优于线性反演方法。但长期以来由于受数学理论发展的限制,人们往往认为地球物理问题的广义非线性反演遵循的规律要比广义线性反演复杂得多,使用迭代方法难以解决问题。本文在非线性最优化理论的最新进展基础上,结合地球物理反演问题,提出广义非线性反演的一系列方法技术。结合层状模型的各向异性弹性参数反演,推导出反演所需的梯度向量,梯度向量中Jacobi矩阵的各个元素的求取利用了与正演完全相同的方法2,并构造一个公式,从单位阵开始来修改函数的海色矩阵的逆矩阵,这样,由于不用求海色矩阵及其逆阵而能够求函数的极小值,减少了计算量。文中给出了实际的反演结果。理论模型和实际资料的反演表明,本方法是求解地球物理反问题的一种有效手段。1 方法原理使用非线性最优化理论导出广义非线性反演方法的第一步是将地球物理反问题转化为一个最优化问题,其最简单的途径是构造如下目标函数f(C)=U(C)-Ud式中:C是由待反演的地下参数构成的参数向量;U(C)为模型响应;Ud是实际的地表观测数据,由C可通过正演求得U(C)。在函数f(C)极小点C3的第k次近似点C(k)处展开f(C)为二次函数f(C)h(C)=f(C(k)+f(C(k)(C-C(k)+122f(C(k)(C-C(k)2式中:f(C(k)为f(C)在C(k)处的梯度;2f(C(k)为f(C)在C(k)处的海色矩阵。以二次函数h(C)逼近f(C),并以h(C)的极小点修正C(k),得到牛顿法的迭代式C(k+1)=C(k)+P(k)其中P(k)=-2f(C(k)-1f(C(k)2f(C(k)-1为f(C)在C(k)处的海色矩阵的逆矩阵。若记g(k)=f(C(k)Gk=2f(C(k)Hk=2f(C(k)-1则有P(k)=-Hkg(k)若f(C)在C(k)处的海色矩阵Gk对称正定,则其逆阵Hk也对称正定,则P(k)为f(C)的下降方向。继续下去有C(k+2)=C(k+1)+P(k+1)其中P(k+1)=-Hk+1g(k+1)而g(k+1)为函数在C(k+1)处的梯度,Hk+1为函数在C(k+1)处的海色阵的逆阵,这样就需用函数f(C)的二阶偏导求Hk+1。广义非线性反演是一个循环迭代、逐步逼近的过程,其主要工作量在计算海色矩阵和求解方程组上。由于地球物理问题一般都较为复杂,计算包含二阶导数的海色矩阵工作量非常大,而海色矩阵是一个与初始猜测模型有关的量,每一次循环迭代过程都必须重新计算,因此减少海色矩阵的计算工作量是提高广义非线性反演速度的关键。本文将给出一个修正Hk的修正矩阵Hk,使Hk+1=Hk+Hk而不必求函数的海色阵及其逆阵。具体构造及证明过程见附录A。2TI介质弹性参数反演2.1 正演问题在反演TI介质弹性参数之前,先研究反演中所用到的正演,本文采用伪谱法作波场正演模拟。在二维TI介质中,地震波波动方程可表示为C6652uy5x2+C5552uy5z2=52uy5t2+fyP52U5X2+Q52U5X5Z+R52U5Z2=52U5t2+?F(1)其中P=C1100C55,Q=0C13+C55C13+C550,R=C4400C33,U=uxuz,?F=fxfz如果令l1=C66525x2+C55525z2-525t2-1l2=P525X2+Q525Z5Z+R525Z2-525t2-1则式(1)可以改写为uy=l1Fy(2a)U=l2F0(2b)式中:Fy和F0为力源函数;L1和L2是非线性算子矩阵。显然,解正演问题的过程就是通过偏微分算子 l1和 l2将力源Fy和F0投影到解析空间的过程。式(2a)表述了正演模拟SH波传播的方程,式(2b)描述了正演模拟P波和SV波耦合在一起传播的方程。212Jacobi矩阵的建立根据常规反演理论可以推知,任何反演方法中,Jacobi矩阵的建立是反演方法成功与否的关键,其中包括Jacobi矩阵建立的基本原理和Jacobi矩阵中各个参数的分布问题,以及Jacobi矩阵中各个参数间的相干性等等问题。根据上面讨论过的TI介质中的弹性波波动方程式(2a)和式(2b),可以得到各个弹性参数相对于各个位移分量之间的偏导数,这些偏导数将组成反演方法中的Jacobi矩阵的元素。根据式(2a)和式(2b),并忽略外力不计,分别用弹性参数作为变量进行偏微分求导,整理后得到下列等式N525x2+L525z2-525t25uy5N=-52uy5x2N525x2+L525z2-525t25uy5L=-52uy5z2N525x2+L525z2-525t25uy5=52uy5t2=1N52uy5x2+L52uy5z2(3)081勘探地球物理进展第30卷P525x2+Q525x5z+R525z2-525t25U5A=-5P5A52U5x2P525x2+Q525x5z+R525z2-525t25U5C=-5R5C52U5z2P525x2+Q525x5z+R525z2-525t25U5F=-5Q5F52U5x5z(4)式中:A,C,F,N和L是用来描述TI介质的弹性参数;为TI介质的密度;uy为在y方向上的位移分量。U=uxuz,5P5A=1000,5R5C=0001,5Q5F=0110如果令r1=-52uy5x2,r2=-52uy5z2,r3=1N52uy5x2+L52uy5z2,R1=-5P5A52U5x2,R2=-5R5C52U5z2,R3=-5Q5F52U5x5z并根据式(1)、式(3)和式(4),则5uy5N=l1r1,5uy5L=l1r2,5uy5=l1r3(5)5U5A=l2R1,5U5C=l2R3,5U5F=l2R3(6)如果用式(5)和式(6)分别与式(2a)和式(2b)相比较,可以发现它们之间有着极大的相似之处,也就是式(2a)与式(5)以及式(2b)与式(6)的表达式基本相同,即它们的解析表达式基本相同,都是通过偏导数差分算子 l1和 l2将震源函数Fy,F0和r1,r2,r3,R1,R2,R3投影到解析空间得到。如果把式(5)和式(6)右边的r1,r2,r3,R1,R2,R3看作是震源项,那么式(5)和式(6)的解就是反演中所需要Jacobi矩阵中的各个参数;并且最为重要的是,利用式(5)和式(6)求解反演中所需要的Ja2cobi矩阵中的各个参数的过程,就是与上述的利用式(2)进行正演过程相同的过程。也就是说,反演中所需求取的Jacobi矩阵中的各个参数,可以通过与正演方法相同的方法求取。它通过对式(5)和式(6)的数值求解即可得到。在求得Jacobi矩阵5U(C(k)5C之后,根据公式g(k)=f(C(k)=2(U(C(k)-Uk)5U(C(k)5C即可求取梯度向量g(k)。利用本文所给的公式(见附录A)修正Hk,其初始矩阵可取单位阵En,即H0En,这样修正Hk+1就不必求函数f(C)的海色矩阵及其逆阵。反演弹性参数的计算程序如下。对于minf(C),xRn,梯度g(k)=f(C(k):取初值C(0),初始矩阵H0=En为单位阵,计算g(0),P(0)=-H0g(0)=-g(0),C(1)=C(0)+P(0),令k=0;求g(k+1),y(k)=g(k+1)-g(k),(Hkg(k+1)Ty(k),若(Hkg(k+1)Ty(k)0时,转步骤;若(Hkg(k+1)Ty(k)0,转,否则转;计 算 Hk=-(Hkg(k+1)(Hkg(k+1)T(Hkg(k+1)Ty(k),H(k+1)=Hk+Hk,P(k+1)=-H(k+1)g(k+1),C(k+2)=C(k+1)+P(k+1);若 g(k+1)0时,需要进行线性搜索,使计算继续),同样可求函数极值。中国地质大学(武汉)能源地球物理研究所招聘启事为加快学科发展与研究所建设,中国地质大学(武汉)能源地球物理研究所现面向社会高薪招聘下列人才(见中国地质大学(武汉)人事处主页)。一、工程岗位工作内容:地震资料数据处理工程师;地震资料地质解释工程师。要求:踏实的工作作风;优良的团队合作精神;8年以上工作经验;处理工程师熟练掌握等地震数据处理系统,解释工程师熟练掌握通用解释软件;负责并承担过大型处理或解释项目。薪资待遇:基本年薪1520万,外加工作量提成(具体面议);该岗位由能源地球物理研究所自行招聘。二、科研与教学岗位学科名称:地球物理学研究方向:地震2测井2地质学;地震学与地震勘探;地震波传播与成像。学科名称:地球探测与信息技术研究方向:地球信息科学与技术。学科名称:构造地质学研究方向:应用构造学。学科名称:能源地质工程研究方向:石油勘探构造分析;沉积学与储层地质学;层序地层学与隐蔽圈闭预测。学科名称:物理学研究方向:波动物理学;光学;声学。学科名称:应用数学研究方向:偏微分方程与反演。学科名称:计算机科学研究方向:计算机软件与理论;软件工程。学科名称:工程力学研究方向:弹性波理论与应用。要求:相关学科博士学位;坚定的科学理念与踏实的工作作风;优良的团队合作精神。薪资待遇:聘任人员享受学校薪酬待遇;另所内待遇面议。应聘者需提供的资料:中、英文简历,近期免冠1寸照片1张;论文发表情况(目录);主要代表作3篇(部);学位证书等其他重要证书材料的复印件。三、联系方式工程岗位应聘者直接电话或E2mail联系。科研与教学岗位应聘者请将应聘材料(书面及电子版)发至下列地址,并请注明申请方向。电话:(027)61063450,(027)62011604电子邮件:wenhuiy 联系人:张老师 於老师地址:湖北省武汉市鲁磨路 中国地质大学(武汉)能源地球物理研究所481勘探地球物理进展第30卷第30卷第3期2007年6月勘探地球物理进展Progress in Exploration GeophysicsVol.30,No.3Jun.,2007ABSTRACTReview of w ave equation prestack depth migration methods.MaShufang,Li Zhenchun.PEG,2006,30(3):153161This paper reviewed the concepts of wave equation pres2tack depth migration and its realizations.Wave equation pres2tack depth migration can be roughly divided into two:single2square2root equation and double2square2root equation.For thesingle2square2root equation migration,the up2going and down2going waves are downward continued respectively,and the cor2relation imaging condition is used.The double2square2root e2quation is based on the survey2sinking concept,the source andthe receiver are downward continued at the same time.Whenthey are on the same place(zero2offset),the wave field valueof zero2time is regarded as the migration value of the space point.True amplitude migration is also discussed,and it is pointed thatthe migration can be realized on common2angle gathers.Theo2retical data processing results showed that wave equation pres2tack depth migration methods are effective tools for complexstructure imaging in media of strong laterally varying velocity.K ey words:wave equation;prestack depth migration;true am2plitude migrationMaShufang,ChinaUniversityofPetroleum,Dongying257061,ChinaApplications of 32D VSP technology and its outlook.Chen Lin.PEG,2006,30(3):16216732D VSP technology can be used to solve the geologicproblems nearby a borehole,and to extract valuable informa2tion about a reservoir.The data obtained by synchronizing 32DVSP and surface seismic acquisition are of great application val2ue,and can be utilizing to optimize survey targets,images thedetailed 32D geology nearby boreholes,and study seismic at2tributes and anisotropic problems.Compared to surface time2lapse seismic,time2lapse 32D VSP is short in acquisition time,cheap,and easy2to2processing,and it is superior in repeatedreservoir characterization,determination of optimal well posi2tion,calculation of reservoir properties in certain areas,moni2toring fluid flow,and reduction of risk in equipment invest2ment.This paper reviewed the applications of 32D VSP tech2nology and outlooked its future trend.K ey words:32D VSP;detailed 32D imaging;reservoir charac2terization;synchronized acquisitionChen Lin,Institute of Geophysical Prospecting,SINOPEC Ex2ploration&Production Research Institute,Nanjing 210014,ChinaEfficacy analysis of harmonic wavelet filtering.Xia Hongrui,G eChuanqing,Dong Jiangwei.PEG,2006,30(3):168174The advantages of harmonic wavelet in detection of vibra2tion malfunction have been accepted,and it has been widelyused in production and scientific research.As an extension ofits reconstruction algorithm,harmonic wavelet filtering has al2so attracted high degree of attention nowadays.In this paper,we introduced the concept of harmonic wavelet and its realiza2tions.Compared to other binary wavelet filtering,harmonicwavelet filtering can overcome the frequency leakage during fil2tering because of its box shape at frequency domain and com2plex filtering.Harmonic wavelet filtering functions same as fi2nite impulse response filtering,and cannot solve the frequencyleakage and fence effect caused by Fourier transform.Some ex2amples were presented in the paper.K ey words:harmonic wavelet transform;binary wavelet trans2form;filtering algorithm;FIR filtering;frequency leakage;fence effectXia Hongrui,Geophysical Prospecting Company of SINOPECJianghan Oilfield Company,Qianjiang 433100,ChinaPhase encoding in prestack migration based on irregular surface.AnQ i,Li Zhenchun,Ye Yueming.PEG,2006,30(3):175178Though wave equation prestack migration produces higherquality images,it has not been widely used in practice due to itslow efficiency.Many attempts addressing this problem focusedon increasing the speed of migration for each individual shot,whereas phase encoding method tries to reduce the number ofshots in migration.Phase encoding method implies horizontalsurface hypothesis.We firstly reconstructed the datum throughwave equation extrapolation,and then implemented phase en2coding in migration at the new datum.The feasibility of theproposed method was tested on theoretical models.K ey words:irregular surface;redatum;phase encoding;pres2tack depth migrationAn Qi,Faculty of Geo2Resources and Information,China Uni2versity of Petroleum,Dongying 257061,ChinaInversion of anisotropic elastic parameters using generalized non2linear inversion method.Sun Yinhang,Yue Youxi.PEG,2006,30(3):179184This paper proposes a generalized nonlinear inversionmethod for anisotropic elastic parameters by introducing non2linear optimization theory into geophysical inversion problems.We derived the gradient vector that is needed in inversion.Theelements of Jacobi matrix in the gradient vector can be calculat2ed through forwardmodeling.We constructed a formulathrough which the inverse of Hessian matrix can be updatedfrom unit matrix.The objective function thus was minimizedwithout knowing the Hessian matrix and its inverse,whichmakes computational cost decreased greatly.Results from realdata inversion were presented.K ey words:anisotropy;elastic parameter;full wave field inver2sion;Hessian matrixSun Yinhang,Faculty of Geo2Resources and Information,Chi2na University of Petroleum,Dongying 257061,ChinaFDTD method for propagation of electromagnetic wave in lossymedia.Li Qingwei,Cao Junxing,Du Xingzhong.PEG,2006,30(3):185189After deduction of difference scheme of Maxwells equa2tion in lossy media,this paper discussed the determination ofparameters in finite difference time domain(FDTD)method indetail.The modifications to the total field boundary conditiondecrease the leakage of incident waves,and lead to no disposalof the four vertices in the total field.This paper discussed theimproved absorbed boundary condition of perfectly matchedlayer(PML),which avoids the splitting of wave components innormal Berenger PML.FDTD forward modeling was carriedout on theoretical models of lossy media.The results depict theuseful information in lossy media,which demonstrates the cor2rectness and reliability of the proposed method.K ey words:finite difference time domain;total field boundarycondition;perfectly matched layer;lossy mediumLi Qingwei,College of Information Engineering,Chengdu Uni2versity of Technology,Chengdu 610059,ChinaStudy of ray multipath and multivalued surface based on elementmodel.Mei Bao,Hu Chaobin,Zhang Shuangjie.PEG,2006,30(3):190193To cope with the strongly folded strata,well developedthrust faults and highly steep structures in the eastern Sichuanarea,this paper put forward a method for 32D geologic modelbuilding based on element model method.The method com2bines 32D geologic modeling with ray theory,and takes full ad2vantage of 32D geologic modeling in describing spatial topologi2cal relations.Methods for dealing with ray multipath and mul2tivalued surface were presented.K ey words:model building;ray trace;element model;data2base;pathMei Bao,Institute of Geophysics and Geomatics,China Uni2versity of Geosciences,Wuhan 430074,ChinaBuilding of geophysical model in piedmont and its application in
