同济大学海洋与地球科学学院博士学位论文地震波数值模拟与反演中几个关键问题研究姓名:董良国申请学位级别:博士专业:固体地球物理学指导教师:马在田20030901
同济大学申请博士学位论文地震波数值模拟与反演中几个关键问题研究董良国摘要地震学的主要任务之一是利用观测到的地震信号提取地下介质物性参数,而叠前地震波形反演比旅行时反演、AVO反演以及叠后反演能获得更多的地下信息。而基于模型的叠前地震波形反演除了包括反演策略和具体的反演方法外,地震波传播理论和地震波传播数值模拟方法在其中占有重要地位。针对叠前地震波形反演中涉及至0的有关问题,本文对其中的一些关键问题进行了研究,其中包括地震波正、反向传播数值计算精度问题、数值计算方法的稳定性问题、数值频散问题、吸收边界问题、起伏地表情况下的自由边界及模拟方法问题、以及反演策略和反演方法问题。地震波正、反向传播波场数值计算精度直接影响着反演的最终结果,而模拟方法和吸收边界条件是决定模拟精度的关键。为了提高地震波正、反向传播数值计算的精度,论文将交错网格数值计算技术和高阶差分方法有机结合,提出了交错网格高阶差分地震波传播数值模拟方法,有效减小了数值频散,提高了地震波模拟的精度,为研究复杂介质模型中地震波传播特征和进行高精度介质参数反演提供了有力工具。另外,论文对这种交错网格高阶差分解法的稳定性进行了理论分析,给出了横向各向同性介质中不同维数、任意时间和空间差分精度情况下统一的稳定性条件,该稳定性条件是对前人关于交错网格差分法地震波模拟稳定性条件论述的高度概括。二维Marmousi模型和三维盐丘模型声波和弹性波模拟实例,充分说明了提出的地震波传播数值模拟方法具有良好的模拟效果。数值频散直接影响着地震波正、反向传播数值计算的精度,因此,论文对有限差分空间和时间离散数值频散这个关键问题进行了理论分析,指出高阶差分法可以在一定程度上降低数值频散、提高计算精度。在吸收边界问题上,论文基于特征波的概念,提出了一种全新的构造声波和弹性波吸收边界条件的特征分析方法,并利用该方法导出了二维和三维声波和弹性波吸收边界条件。和其它方法相比,该方法的主要特点是可以方便地构造各边角区域的吸收边界条件,并和各面的吸收边界条件具有相同的吸收精度,保证了地震波正、反向传播波场计算的稳定,并很好地保持了边界区域地震波形的一致性。起伏地表地震波传播数值模拟是波场模拟中的一个难点,适合地形起
摘要伏的模拟方法以及正确的自由边界条件实现方式是其中的关键问题。为此,论文从自由边界的概念出发,得到了地表起伏情况下用应力表示的自由边界条件,利用纵向坐标变换思路将曲网格变换为矩形网格,为利用提出的交错网格高阶差分数值模拟方法铺平了道路。论文同时将第四章提出的吸收边界条件推广到起伏地表。地下介质物性参数提取是地震学的重要目标,鉴于地震反演强烈的非线性和多解性问题,论文在对反射地震波最小平方反演目标函数性态进行分析的基础上,提出了分步骤、多尺度反演策略。利用深度偏移结果构造目标函数,通过遗传算法反演速度参数变化的低频趋势,以此为初始速度模型,利用共轭梯度叠前波形反演求取速度空间变化的高波数成分。论文在梯度方向的计算中,通过引进不同炮点、不同检波点的权函数,提高了梯度方向的计算精度,加快了反演的收敛速度。采用地震数据频率多尺度反演策略,将上述地震波传播高阶差分模拟方法应用到参数反演之中,提高了反传播类地震波形反演的精度,并对地震数据信噪比、初始模型、地震子波对反演结果的影响进行了数值分析。通过将上述反演方法和思路应用于南海实际资料中,得到了海底精细一维速度结构,。上述反演技术如果要大量应用于实际,还需要在子波提取等方面进一步改进,但2D和3D声波、弹性波数值模拟技术已形成软件,并在全国许多油田和地球物理勘探单位得到广泛应用。Ⅱ
ABSTRACTOneoftheobjectivesofseismologyistoobtaintheea/'thpropertiesusingtheobservedseismicdata.Moregeologicalinformationcanbeobtainedbypre..stackseismicwaveforminversionthanbytravel-timeinversion,AVOinversionandpost.stackwaveforminversion.Apartfromtheinversionstrategyandmethods,theseismicwavepropagationtheoryandsimulationtechniquesplayanimportantroleinthemodel·basedpre·stackseismicwaveforminversion.AccordingtOtheaspectsintheseismicwaveforminversion,severalkeyproblemsareanalyzedandstudied,includingthemodelingmethodsandtechniquesofseismicwavepropagation,theaccuracyandstabilityofthenumericalmethods,thefreesurfaceboundaryconditionsandsimulationmethodinmodelswithtopography,theabsorbingboundaryconditions,andwaveforminversionstrategiesandtechniques.Theaccuracyofseismicwaveforwardandreversepropagationcalculationseriouslyaffecttheinversionresults.Modelingmethodsandabsorbingboundaryconditionsaretwokeypointsinthenumericalsimulationofseismicwavepropagation.Inordertoincreasetheaccuracyanddecreasethenumericaldispersion,thestaggered·grid,high-orderdifferencemethodisproposedtosimulatetheseismicwavepropagation,whichcombinedthestaggered一鲥dtechniqueandthehip.orderdifferencemethod.ThismethodCanbeusedinthestudyofseismicWaVepropagationandinversion。Furthermore,thestabilityofthismethodiSstudied.AuniquestabilityconditioniSderivedfordifferentdimensions,anddifferenttemporalandspatialdifferenceordersintransverselyisotropicmedia.Itisasummaryofstabilityconditionstothisstaggered—gridmethod,proposedbypreviousauthors.Theexamplesof2-DMarmousimodeland3-DSaltmodelshowthatthismethodCallproducefairlygoodsimulatedresults.Numericaldispersiondecreasestheaccuracyofseismicwaveforwardandreversepropagationcalculation.Thefinite—differencespatialandtemporaldispersionsaretheoreticallyanalyzed,andIdrawtheconclusionthatlli曲.orderfinite.differencecalldecreasethespatialandtemporalgriddispersion.Anewapproachtoconstructabsorbingboundaryconditions(ABCs)foracousticandelasticwaveequations(forbothisotropicandakindofanisotropicmedia(VTI))m
ABSTRACTispresented.Theabsorbingboundaryconditionsforbothacousticandelasticwaveequations(2-Dand3一D)arcalsoderived.,nlisappmachusesthecharacteristicanalysismethodtoseparateseismicwavesintooutgoingandincomingwaves.ApartfromtheABCsatthefourlinesin2一dimensionandsixsurfacesin3-dimension,theABCsatfourcornersin2-dimensionandeightcomerSand12linesin3-dimensionareconvenientlyconstructed.DifferentboundaryregionsaredealtwithrespectivelySOthatthestabilityofthenumericalalgorithmcailbeachievedintheseismicwavesimulation.Anannoyingproblemintheseismicwavemodelingisthesurfacetopography.Appropriatesimulationmethodandgoodimplementationoffree-surfaceconditionarcrequiredtodealwimthisproblem.Anexplicitfleesurfaceboundaryconditionfortopographyisderived.Usingtheideaofcoordinatetransformation,thecurved鲥distransformedtorectangulargridSO’thattheproposedstaggered-gridhigh-orderdifferencemethodCanbeusedintheseismicwavesimulationinthepresenceoftopography.Meanwhile,theabsorbillgboundaryconditions,proposedinchapter4,aleextendedtothesituationofsurfacetopography.Owingtothenonlinearpropertyoftheseismicinverseproblem,thecharacterofthefitnessoftheleastsquaresinverseproblemisexamined,andthenthemulti-stepandmulti—scaleinversionstrategyispresented.TheGAalgorithmisusedtoinvertthelongwavelengthpartofthevelocity,withthedepthmigratedresultstoformthefitness.Thentheperturbationsofthevelocityareobtainedusingpre-stackseismicwaveforminversion,wi也theCGiterativealgorithm弱thet001.Inthecalculationofthegradientdirection,differentweightsfordifferentshotsanddifferentreceiverpositionsare眇cdt0raisetheaccuracyandincreasetheinverseefficiency.Multi-frequencybandstrategyandtheproposedstaggered—gridhigh-orderdifferencesimulationmethodaleusedtoraisetheinverseresults.Howthes/nratioofseismicdata,initialmodelandseismicwaveletaffecttheinverseresultsiSalsohumericaltested.Finally,thefinevelocitystructure(1-D)ofasiteinSouthChinaSeaisobtainedusingthisinverseprocedure.IV
同济大学中请博士学位论文地震波数值模拟与反演中几个关键问题研究董良国第一章前言1.1地震波正反演研究内容地震学的根本任务是根据观测到的地震信号研究和提取有关地下介质的物性信息,如速度、密度等。这一问题的解决有正演和反演两种途径(vijayDimri,1992),它们也是弹性动力学研究的两个基本方面。正问题是研究在给定震源和介质性质时地震波的传播规律,而反问题则是研究如何根据观测到的地震信号反演地下介质的有关性质。对天然地震学而言,反问题还包含反演产生波场的震源特征。地震波传播正演是我们认识地震波传播规律的重要途径,是利用地震波研究地球内部结构、勘探地下自然资源、进行地下工程探测的理论基础。在弹性理论的基础上,在线性小应变的前提下,地震学正问题归结为初值条件的线性偏微分方程的定解问题。由于一般情况下无法获得该类定解问题的解析解,所以如何获得这类偏微分方程的数值解,即如何进行地震波传播数值模拟成为正问题研究的一个核心问题。由于在医学CT、无损检测、地震预报、能源勘探、海洋声学、探地雷达以及利用观测到的地震数据确定地球内部结构等方面的诸多问题最终都归结为弹性动力学反问题,因此地震学中弹性动力学反问题的研究具有重要的理论意义。在油气勘探方面,目前我国的岩性油气藏所占比重越来越大,目前已经超过了50%,而常规成像分辨率达不到识别岩性体的要求,这些岩性体对地震波的影响主要表现在波形上,决定了地震反演在油气勘探中的作用越来越大,因此,地震学中反问题的研究还有广泛的应用价值。正因为如此,弹性动力学反问题的研究,包括反演理论研究以及应用研究,在近几十年来一直是数学家和地球物理学家研究的一个热点。由于该问题的非线性、不适定性和解的非唯一性,加上地震数据的带限性、观测区域的有限性以及噪音的干扰,造成反演问题的求解非常困难。地震学反问题一般可以归结为非线性泛函的极小问题(刘家琦,1983),这一问题的解决主要决定于三个方面:地震波传播理论、地震数据观测技术和地震波反演方法。从历史上看,野外观测技术的每次改进(如数字地震仪、多次覆盖技术、三维观测技术),都会引起勘探技术革命性的飞跃,造成探测精度大幅度提高。相信形象再现地震波传播特征的数值模拟技术的发展,也会有效提高反演的精度。而要根本解决地下介质参数问题,反演理论还有待突破。可以预见,地震学反问题重大突破要依赖于地震波传播理论与非线性理论的有机
第一幸前言结合(杨文采,1995,2002)。其中,优化过程中地震波正演占据举足轻重的地位,地震波传播正演精度一定程度上决定了反演的精度,而且反演效率主要取决于正演方法的效率,因为在目前的地震波形反演方法中,地震波正演占据了整个反演计算量的95%以上!因此,地震波反演离不开正演。本论文对地震波传播数值模拟方法以及叠前地震波形反演方法进行了研究,重点研究了其中的几个关键技术问题。1.2地震波模拟与反演研究现状1.2.1地震波模拟研究现状地震波正演贯穿于地震学发展的始终,其中预测给定地质模型中特定点上地震响应的地震波传播模拟是正演的一个重要方面。物理模拟和数值模拟是地震波传播模拟的两种不同方式。物理模拟是在实验室内,根据一定的相似比制作物理模型,激发的超声波在模型中传播,通过记录这些超声波来研究波在特定介质中的传播规律。整个模拟过程是客观的物理过程,但目前模型制作工艺跟不上我们对复杂介质描述的要求,换能器性能上也存在缺陷(如低频换能器制作困难、指向性差等),模拟工期较长,模型修改困难,费用较高,这些缺点严重阻碍了物理模拟技术在实际中的广泛应用。地震波传播数值模拟是在一定的地震波传播理论基础上,通过数值计算来模拟地震波的传播。尽管数值模拟技术所依据的地震波传播理论是实际地震波传播规律的不同程度的近似,但该技术灵活方便、模拟快速等优点使其在理论研究和实际中得到广泛应用,近年来计算机技术的日新月异进一步推动了数值模拟技术的发展。所以,地震波传播数值模拟是地震数据处理的基石(Carcione,2002),也是研究地震波传播理论、指导地震数据采集、进行地震资料解释的有力工具,同时又是许多地震反演算法中的一个关键步骤。地震波数值模拟方法是和地震波传播理论紧密结合在一起的。由于可以用射线理论、积分方程、微分方程来描述地震波的传播,模拟方法也相应地有射线追踪方法、积分方程数值求解方法以及微分方程数值求解方法。射线追踪方法通过求解程函方程计算地震波旅行时,通过求解传输方程计算地震波振幅。该方法以高频近似为前提,适合于物性缓变模型中地震波传播模拟。模型简单时该方法具有计算速度快的突出优点,正因为如此,它在地震成像、旅行时层析等方面得到广泛应用。也正是由于其高频近似,该方法不适合物性参数变化较大模型中地震波的传播模拟。2
同济大学申请博士学位论文地震波数值模拟与反演中几个关键问题研究董良囤积分方程数值求解地震波数值模拟方法是基于惠更斯原理而得到的一种波‘场计算方法,它又可以分为体积分方法和边界积分方法。由于该方法的半解析特征,使其在成像、反演理论研究和公式推导方面具有得天独厚的优势(有关这方面的应用可参见Bleistein(1987)的专著)。由于涉及Green函数的计算,该方法一般适合于模拟具有特定边界地质体产生的地震波,而要求该地质体周围为均匀介质。因此,该方法的适应范围受到严格限制。微分方程方法是对计算区域网格化,通过数值求解描述地震波传播的微分方程来模拟波的传播。就目前看来,该方法对模型没有任何限制,在地震波模拟中使用最为广泛,主要问题是计算量比较大,对计算机内存要求较高。其中,有限差分法(FD)、有限元法(FE)以及傅立叶变换法(PS)是这类模拟方法中使用较多的方法。近年来还出现的界于有限差分法和有限元方法之间的有限体方法(FV),在理论上应该具有有限元方法网格剖分的灵活性,又具有有限差分法计算快速的特点,但在简单的矩形网格情况下,该方法完全退化为有限差分法。上述方法具有各自的优缺点,有限元法(邵秀民,蓝志凌,1995)剖分灵活,适合于处理不同形状的介质分界面以及起伏的自由界面,但使用不方便(主要是灵活剖分的网格点上的物性参数难以方便设定),效率较低,所需内存较大。傅立叶变换法(KosloffandBaysal,1982;Reshefeta1.,1988)主要思想是利用Fourier变换计算波场的空间导数,计算精度高,但不适合复杂模型中地震波模拟,边界处理也比较困难。与Fourier变换法思路相似,GottliebandOrszag(1977)提出的一种谱方法利用Chebyshev多项式拟合来计算波场的空间导数,提高了模拟的精度,但复杂模型的变网格描述不方便,限制了它在实际中的应用。有限差分法(Kellyeta1.,1976)使用方便灵活,应用比较广,目前国际勘探地球物理界著名的二维Marmousi模型以及三维SALT模型的地震数据合成均是利用FD方法。在2002年SEG年会上推出的MarmousiII模型也是利用差分法完成的(GaryMartin,2002)。自从Altermaneta1.(1968)将FD法引入到地震波模拟中以来,如何进一步提高模拟精度,减小数值频散,一直都是有限差分法研究的重点。从二阶差分(Alford,eta1.,1976)到高阶差分(Dablain,1986),从规则网格到交错网格(Virieux,1986;OzdenvaranfMcMechan,1997),其目的都是这样。然而,在许多情况下,这一问题并没有得到根本解决,尤其是模拟低速介质中的弹性波传播时,数值频散问题仍很严重。另外,模拟区域的边界处理也一直是差分法及其它网格模拟方法必须面对的问题,这其中包括介质分界面以及地表自由边界。尽管发展的多种吸收边界条件对人为边界反射都具有不同程度的吸收能力,但没有一种方法能够将人为边界反射完全吸收,这些反射回模拟区域的人3