二维频率域声波方程正演模拟

合集下载

二维常密度声波波动方程

二维常密度声波波动方程
∂^2p/∂t^2 = c^2(∂^2p/∂x^2 + ∂^2p/∂y^2)。

在这个方程中，p代表声压，t代表时间，c代表声速，x和y 代表空间中的两个方向。

这个方程描述了声波在二维空间中随时间和空间的演变。

声波在传播过程中会受到空气密度、压力等因素的影响，而这些因素会影响声波在空间中的传播速度和方式。

从物理角度来看，这个方程可以用来描述声波在二维空间中的传播规律，可以帮助我们理解声波的传播特性和行为。

从工程角度来看，这个方程可以应用于声学领域的相关工程问题，比如声波传感器设计、声场模拟等方面。

从数学角度来看，这个方程是一个偏微分方程，可以通过数值方法或解析方法进行求解，因此在数学建模和分析中具有重要意义。

总的来说，二维常密度声波波动方程是描述二维空间中声波传播规律的重要方程，对于理解声波的物理特性、工程应用以及数学分析都具有重要意义。

二维弹性波方程频率空间域有限差分正演研究

第１Ｏ卷第６期
２Ｏ１３年１１月
工往球物理号旅
ＣＨＩＮＥＳＥＪ０ＵＲＮＡＬＯＦＥＮＧＩＮＥＥＲＩＮＧＧＥＯＰＨＹＳＩＣＳ
Ｖ０ＬｌＯ，Ｎｏ．６
ＮＯＶ．，２０１３
文章编号：１６７２ —７９４０（２０１３）０６一Ｏ８５８一Ｏ５
ＮａｎｃｈａｎｇＪｉａｎｇｘｉ３３０２０１，Ｃｈｉｎａ）
Ａｂｓｔｒａｃｔ：Ｔｈｅｔｉｍｅ — ｓｐａｃｅｄｏｍａｉｎｆｉｎｉｔｅｄｉｆｆｅｒｅｎｃｅｆｏｒｗａｒｄｍｏｄｅｌｉｎｇｉｓｒｅｌａｔｉｖｅｌｙｍａｔｕｒｅ，ｂｕｔｔｈｉｓｍｅｔｈｏｄｈａｓｔｒｕｎｃａｔｉｏｎｅｒｒｏｒｓｗｈｉｃｈｃｏｎｓｔａｎｔｌｙａｃｃｕｍｕｌａｔｅｗｉｔｈｔｈｅｉｎｃｒｅａｓｅｄａ — ｍｏｕｎｔｏｆｃｏｍｐｕｔａｔｉｏｎ．Ｈｏｗｅｖｅｒ，ｉｎａｃｃｏｒｄａｎｃｅｗｉｔｈｔｈｅｆｒｅｑｕｅｎｃｙｓｌｉｃｅｓ，ｆｒｅｑｕｅｎｃｙ—
３．江西省地质矿产勘查开发局物探化探大队，江西南昌３３０２０１）
摘要：时间空间域有限差分正演模拟方法比较成熟，但是，该方法存在截断误差，并且该误差会随着计算

CSC频率—空间域波动方程数值模拟

CSC频率—空间域波动方程数值模拟吕晓春;顾汉明;成景旺;周丽【摘要】针对频率空间域波动方程数值模拟需要巨大内存空间的现状,提出了利用列索引压缩存储(CSC)技术存储大型稀疏非对称复数型的矩阵系数.CSC技术将系数矩阵转化为三个一维数组来存储,分别存储系数非零元素、非零元素对应所在的行以及每列起始非零元素所在位置.经CSC技术压缩存储后显著减少了内存空间及计算量,在计算时只有少许的非零元素参加计算,且根据三个一维数组可以简便地找到对应的非零元素,进而采用LU分解快速而精确地求解.本文基于Jo等提出的最优化9点差分方法,首次应用CSC技术在频率空间域进行二维声波方程数值模拟.通过对Corner-edge模型和二维Marmousi模型进行试算,可以显著节省内存需求,明显提高计算速度,进而得到精度较高的正演结果.【期刊名称】《石油地球物理勘探》【年(卷),期】2014(049)002【总页数】7页(P288-294)【关键词】频率—空间域;CSC;系数矩阵;一维数组;LU分解;数值模拟【作者】吕晓春;顾汉明;成景旺;周丽【作者单位】中国地质大学(武汉)地球物理与空间信息学院,湖北武汉430074;中国地质大学(武汉)构造与油气资源教育部重点实验室,湖北武汉430074;中国地质大学(武汉)地球物理与空间信息学院,湖北武汉430074;中国地质大学(武汉)构造与油气资源教育部重点实验室,湖北武汉430074;中国地质大学(武汉)地球物理与空间信息学院,湖北武汉430074;中国地质大学(武汉)构造与油气资源教育部重点实验室,湖北武汉430074;中国地质大学(武汉)地球物理与空间信息学院,湖北武汉430074【正文语种】中文【中图分类】P6311 引言为了反演地下介质参数或研究地震波在各种复杂介质中的传播机制，需要进行波场数值模拟。

波动方程的数值模拟方法包括有限差分法、有限元法、伪谱法等。

在时间—空间域的数值模拟技术已较成熟，并广泛应用于复杂介质的正演模拟中。

二维地震正演模拟方法技术研究

二维地震正演模拟方法技术研究一、本文概述随着地球物理学的深入发展和油气勘探的不断推进，二维地震正演模拟方法技术在地震勘探领域的应用越来越广泛。

该技术通过模拟地震波在地下介质中的传播过程，为地震资料解释、储层预测和油气勘探提供重要的理论支撑和实践指导。

本文旨在深入研究二维地震正演模拟方法技术，探讨其基本原理、发展历程以及当前的研究热点和难点，为进一步提高地震勘探的精度和效率提供理论支持和技术保障。

本文将对二维地震正演模拟方法技术的基本概念进行阐述，包括其定义、特点以及应用领域等。

接着，回顾二维地震正演模拟方法技术的发展历程，分析其在不同阶段的主要特点和优缺点。

在此基础上，重点探讨当前二维地震正演模拟方法技术面临的主要挑战和难点，如复杂地质条件下的模拟精度问题、大规模计算的效率问题等。

针对这些挑战和难点，本文将进一步分析现有的解决方案和发展趋势，如基于高性能计算的并行化技术、基于人工智能的反演方法等。

同时，结合具体的应用案例，分析二维地震正演模拟方法技术在油气勘探、矿产资源调查等领域的实际应用效果，以验证其有效性和可靠性。

本文将对二维地震正演模拟方法技术的未来发展进行展望，提出可能的研究方向和应用前景。

通过本文的研究，旨在为推动二维地震正演模拟方法技术的发展和应用提供有益的参考和借鉴。

二、二维地震正演模拟理论基础二维地震正演模拟是地球物理学中一种重要的方法，其理论基础主要基于波动方程和地震波的传播原理。

在二维空间中，地震波的传播受到介质速度、密度、弹性等因素的影响，这些因素决定了波场的空间分布和时间变化。

理解和应用波动方程是二维地震正演模拟的关键。

波动方程是描述波在介质中传播的基本方程，对于地震波而言，常用的波动方程有弹性波方程和声波方程。

在二维正演模拟中，我们通常采用声波方程，因为它相对简单且能够较好地模拟地震波的主要特征。

声波方程描述了声波在弹性介质中的传播规律，包括波速、振幅、相位等参数的变化。

二维声场仿真实验报告

二维声场仿真实验报告实验目的：本实验旨在通过二维声场仿真，研究声音在不同环境中的传播规律，并探索不同因素对声音传播的影响。

实验装置及原理：本实验使用MATLAB软件进行二维声场仿真。

首先，确定二维空间中声源的位置、声源频率、声源强度等参数。

其次，根据环境的声学特性设置相应的边界条件和传播介质特性。

然后，通过声学波动方程的数值解法，计算出声场中各点的声压值。

最后，根据计算结果绘制声场图像，以便观察声音在空间中的传播规律。

实验步骤：1. 确定声源参数：包括声源位置、声源频率、声源强度等，根据实验要求进行设置。

2. 设置边界条件：根据实验环境的特性，设置声场的边界条件。

例如，可以选择开放空间、封闭空间或者带有障碍物的空间等。

3. 设定传播介质：确定声场的传播介质，包括密度、声速等参数。

根据实际情况进行设置。

4. 利用声学波动方程进行仿真计算：根据声场的边界条件和传播介质特性，利用声学波动方程的数值解法进行仿真计算，计算出声场中各个点的声压值。

5. 绘制声场图像：根据仿真计算结果，绘制声场的声压等值线图或3D图像，以便观察声音在空间中的传播规律。

实验结果与讨论：根据实验设定的参数，进行了二维声场仿真实验，并获得了声场的声压分布图像。

通过观察图像，可以发现声音在不同环境中的传播规律。

例如，在开放空间中，声音的传播相对均匀，声压等值线比较均衡；而在封闭空间中，声音受到边界的反射和干扰，声压等值线出现了不规则的分布。

在实验过程中，还对不同因素对声音传播的影响进行了讨论。

在传播介质方面，不同密度和声速的介质会对声音的传播速度和路径产生影响。

在边界条件方面，开放空间和封闭空间的边界反射特性不同，会导致声音传播效果的差异。

结论：通过本实验的二维声场仿真，我们可以更好地理解声音在空间中的传播规律，并了解不同因素对声音传播的影响。

这对于优化声音传播环境设计、改善音响系统效果等具有一定的指导意义。

同时，通过MATLAB软件进行声场仿真，也展示了数字仿真技术在声学研究中的应用前景。

声波波动方程正演模拟分析研究

Δt
c
在数值计算中,生成的强边界反射会对中心波
其中:
σ0 =l
og
1
R
3
x
δ
2
(
11)
3Vp
,
R 是理论反射系数;
δ
2δ
是 PML 的厚度;Vp 为速度;在此基础上可推演由
PML 边界条件进行交错网格的有限差分格式。
理论上,
PML 法对各种入射角和频率下的地震
弥散,产生数值频散现象 [9],严重影响正演模拟的精
度。为了减轻数值频散,通常可采用以下方式:① 调
整恰当的时间和空间离散步长,尤其空间步长不宜
过大过小,而时间步长相对越小越好,但会受到实际
计算效率的限制。 ② 通过提高差分阶数,其中提高
空间差分阶数易实现且能有效降低频散,但随着空
2023 年 7 月
第 13 期总第 527 期
Ju
l
y2023
No.
13 To
t
a
lNo.
527
内蒙古科技与经济
I
nne
r Mongo
l
i
aSc
i
enc
eTe
chno
l
ogy & Ec
onomy
声波波动方程正演模拟分析研究
朱晓洁
(中国石化胜利油田分公司海洋采油厂,山东东营 257237)
波的吸收效果良好,吸收效率强于传统的吸收边界
条件法 [14]。因此在计算区域加入吸收边界后,可以

波动方程正演在地震勘探设计中的应用

一
① 中国教育；ｔｌｔ［Ｎ］．２００７ — ０９。
参考文献：
［１】劳动社会保障部职业技能鉴定中心组《职业社会能力》训练手册【Ｍ】．北京；人民出版社，２００８．［２】解志斌．高职思想政治理论课职业核心能力项目化教改的实证考量Ｕ］．牡丹江大学学报，２０１０，（８）．［３】李飞．浅谈思想政治理论课教学模式与高职学生核心能力培养Ｕ１．中国轻工教育，２０１１，（６）．
越重要，利用现有的地震资料，结合地震反演建立合理的地质模型，运用波动方程正演技术设计地震勘探采集方法．研究偏移成像效果，分析采集结果能否达到地质目标的需求，通过反复模拟处理分析，优化观测方案，使得最终确定的采
集方法能够解决地质问题，获得高质量的地震资料，减低勘探风险。
采用高阶差分法得到声波方程的数值解为：
ｃ筹：去ＩＶＩ
。
（ △ｘ２Ｍ）
高阶差分系数：
１２３
４４４
…
Ｍ。
４
ｃ【Ｍ）
（Ｍ）
二、波动方程原理
ｌ２３
６６６
… Ｍｃ２
一
Ｏｘ
Ｏｘ
）＋茜（吉等）＋ｏ（１０ｚ）＝亡鲁
（ｘ＋ｍ△ｘ）＿２ｆ（ｘ）＋ｆ（ｘ — ｍ△ ｘ）］＋０
（ｘ，Ｙ，ｚ，ｔ）８（Ｘ－Ｘｏ，Ｙ－Ｙｏ，Ｚ－Ｚ０）

频率多尺度全波形速度反演

频率多尺度全波形速度反演张文生;罗嘉;滕吉文【摘要】以二维声波方程为模型,在时间域深入研究了全波形速度反演,全波形反演要解一个非线性的最小二乘问题,是一个极小化模拟数据与已知数据之间残量的过程.针对全波形反演易陷入局部极值的困难,本文提出了基于不同尺度的频率数据的“逐级反演”策略,即先基于低频尺度的波场信息进行反演,得出一个合理的初始模型,然后再利用其他不同尺度频率的波场进行反演,并且用前一尺度的迭代反演结果作为下一尺度反演的初始模型,这样逐级进行反演.文中详细阐述和推导了理论方法及公式,包括有限差分正演模拟、速度模型修正、梯度计算和算法描述,并以Marmousi复杂构造模型为例,进行了MPI并行全波形反演数值计算,得到了较好的反演结果,验证了方法的有效性和稳健性.【期刊名称】《地球物理学报》【年(卷),期】2015(058)001【总页数】13页(P216-228)【关键词】声波方程;频率多尺度;时间域;全波形反演;速度;BFGS;Marmousi模型;有限差分法;MPI并行【作者】张文生;罗嘉;滕吉文【作者单位】中国科学院数学与系统科学研究院,计算数学与科学工程计算研究所,LSEC,北京100190;中国科学院数学与系统科学研究院,计算数学与科学工程计算研究所,LSEC,北京100190;中国科学院地质与地球物理研究所,北京100029【正文语种】中文【中图分类】P3151 引言地震全波形反演是利用地表或钻孔中观测到的叠前波场记录信息，推测地球内部介质的物性参数，以确定含油气构造.速度是重要的物性参数之一，速度反演或建模是地球物理偏移成像和地震资料解释的重要基础.走时反演或成析成像（例如，Fawcett et al.，1985；Dahlen et al.，2000；Hung et al.，2004；Gautier et al.，2008；Tong et al.，2011）是速度反演的一种重要方法，但该方法只利用了波场的到时信息，本质上是一种高频近似方法，适于对区域背景速度场成像，对复杂速度构造模型成像精度不足.全波形反演（例如，Tarantola，1984；Mora，1987；Yang，1993；Pratt et al.，1998）利用波场的走时、振幅和相位信息来反演，具有更好的精度和更高的分辨率，但计算量大，难度大，特别是二维和三维情形.从20世纪80年代以来，人们一直致力于这方面的研究，发展了多种地震全波形反演方法，这些方法相互交叉或综合，并无严格界限，大致可如下分类.基于求解波动方程计算域的不同，可分时间域反演（例如，Tarantola，1984；Bleistein，1987；Zhang et al.，2013）、频率域反演（Pratt et al.，1988；Song，1995；Pratt et al.，1998；Prat，1999；Sirgue et al.，2004）和Laplace域反演（Shin et al.，2008，2009；Ha et al.，2012）.基于求解方式的不同，有局部线性化梯度类迭代反演（Tarantola，1984，1987；Gauthier et al.，1986；Bleistein，1987；Mora，1987）、混沌反演（Yang，1993）和全局搜索优化反演等.线性反演不能恢复速度模型的长波场分量（Devancy，1984），波动方程反问题的主要方法是非线性优化迭代反演方法和Born反演（Bleistein，1987）.Tarantola首先（Tarantola，1984，1986）对完全非线性全波形反演进行了详细的研究，数学上归结为一个局部最优化问题.求解这类问题时，扰动模型沿目标函数的梯度搜索，其中梯度可由震源激发的入射波场和残差反传播后的波场计算得到，这种方法可以避免直接通过模型参数扰动来计算Frechet导数，大大减少了计算量，使得二维时间域波动方程全波形反演成为可能，该思想已用到声波和弹性波方程的全波形反演中（Gauthier et al.，1986；Mora，1987）；此外，Sheen等（Sheen et al.，2006）对弹性介质在时间域用高斯牛顿法进行了全波形反演，Choi等（Choi et al.，2008）对声学弹性复合介质中的多分量数据用共轭梯度法进行了二维全波形反演，Epanomeritakis等（Epanomeritakis et al.，2008）结合高斯牛顿法和共轭梯度法对三维弹性波方程进行了全波形反演.将时间域的声波方程和弹性波方程变换到频率域，然后对给定的频率再进行正演模拟，反演也在频率域中求解（Pratt et al.，1988；Song，1995；Pratt et al.，1998；Prat，1999；Sirgue et al.，2004），这就是频率域全波形反演.Pratt等（Pratt et al.，1988）首先将全波形反演思想用于频率域中，频率域反演和时间域反演的思想是完全一致的，只是频率域的波场是互相解耦的，可以对每个频率单独进行反演，频率域反演对初始模型也存在较大的依赖性.类似地，全波形反演还可在Laplace域中进行（Shin et al.，2008；Shin et al.，2009；Ha，2012）.全波形反演数值上常用迭代法来求解，总体上是一个梯度类的优化迭代求解过程，有精度高的优点，但对初始模型有严重依赖性，这是梯度类反演方法的共性.Gauthier等（1986）曾证明当初始速度模型与真实模型有10%的误差时，全波形反演就失效.在这种情况下，梯度类方法趋向找到局部极值点.为了提高全波形反演对初值的稳健性，Bunks等（1995）提出了多尺度全波形反演方法，通过将问题分解为不同尺度求解以保证反演过程稳定收敛.此外，人们还研究了用全局优化搜索的方法来求解.全局优化搜索基于随机过程，不需目标函数的导数信息，也不考虑目标函数的形状，不需要一个好的初值模型，而且不会陷入局部极值点，例如，Jin等（1994）提出了蒙特卡洛法反演，Smabridge等（1992）提出了遗传算法反演，Sen和Stoffa（Stoffa et al.，1991；Sen et al.，1995）提出了模拟退火法.全局优化搜索由于为求得目标函数的全局极值所导致的计算量极大，特别是二维和三维反演.依据线性化反演理论，当一个初始模型在观测数据和模拟数据之间引起大的动力学残差时，速度模型的扰动对残差的平方就没有影响，从而数据对模型参数的梯度为零，但这时求得的极值点显然不是全局极值点，这就造成速度反演的大量极值点，不但难以求得好的收敛结果或甚至发散，而且反复迭代计算也导致巨大的计算量.为此，本文研究了时间域频率多尺度全波形反演方法，提出了基于不同尺度的频率数据的逐级反演策略，取得了明显效果.2 理论方法2.1 有限差分正演模拟在时间域上，考虑二维声波方程：其中x为地面横向坐标，z为深度纵向坐标，u（x，z，t）表示压力，v（x，z）为介质速度，f（t）是震源函数，（xs，zs）为震源位置.若介质在激发前静止，则初始条件为有限差分法是数值模拟波场传播的重要方法之一，具有计算效率高的优点.本文用有限差分法来求解波动方程.设Δx和Δz分别为x和z方向的离散步长，Δt为时间步长，用表示lΔt时刻及空间位置为（nΔx，mΔz）处的波场值，vn，m 表示空间位置（nΔx，mΔz）处的速度离散值.用二阶中心差分格式对方程（1）进行离散，其中，（ns，ms）为震源坐标，这里记M＝NxNz为空间离散点数，在反演波速时也即优化问题的未知量个数.初始条件（2）的离散格式为应用离散Fourier级数法分析易得，该二阶差分格式的数值稳定条件为由于计算都是在有限区域上进行的，为了消除边界反射模拟无界区域中的波场传播，需要应用吸收边界条件（Clayton et al.，1977；Berenger，1994；Yang et al.，2002，2003）.通常计算区域都是规则的矩形域，不妨设为Ω＝（0，X）×（0，Z），我们采用局部化的二阶边界吸收条件（Clayton et al.，1977）：相应的差分格式是：其中，显然，边界条件的差分格式（11）—（14）也是二阶精度.2.2 模型修正全波形反演是一个迭代极小化观测数据与拟合数据之间的残量即目标函数来反演速度v的过程.残量使用l2范数来度量，目标函数可写成：其中，ucal表示数值模拟的波场值，uobs表示已知的地震记录数据.模型参数v 表示整个离散区域的速度，是一个M维向量.目标函数J（v）的极小化是从一个给定的初始值v0出发进行搜索，这是一个局部优化问题.由于数据与模型参数间的非线性关系，需要多次迭代才能收敛到目标函数在v0附近的局部极小点.每次迭代都寻找一个下降方向p和搜索步长α，以确保新的迭代步满足即满足目标函数值下降，注意这里p也是一个M维向量.步长α可用一维线搜索方法（袁亚湘等，2005）得到.计算搜索方向的方法有很多，一般都要基于求解目标函数的一阶导数信息，即梯度；但单纯的梯度法在解附近收敛缓慢，需要用目标函数的二阶导数信息，即Hessian矩阵来修正搜索方向，这类通称为牛顿类方法.为此，将目标函数J（vk＋p）在迭代点vk处作泰勒近似展开，忽略高阶项，得对上式关于分量pl求偏导，得目标函数当时达到极小值，所以有方向p的表达式实际计算中用拟牛顿法来计算，如DFP法（Fletcher et al.，1963；Davidon，1991）和BFGS法（Broyden，1970；Fletcher，1970；Goldfarb，1970；Shanno，1970），即不直接计算式（23）中的二阶导数，而是在每个迭代步用某个对称非奇异且易求逆的矩阵Bk代替该二阶导数矩阵.实际计算表明，BFGS方法是拟牛顿法中最有效的一个方法，Bk和其逆均有递推的修正公式.对传统的BFGS方法，由于Hessian阵初值选取不当或者当前迭代步Hessian阵条件数很坏，会出现收敛慢的情况.为了克服这种情况，可在更新Hessian阵之前对矩阵进行一个尺度变换，即用τkBk代替Bk，τk为变尺度因子，这种变尺度的方法能加速BFGS法在当前迭代步的下降率（Nocedal et al.1993）.在具体计算中，为节省计算内存，常采用有限内存的BFGS方法即L－BFGS方法来计算.下面给出全波形反演的算法流程.步骤1：给定初始模型v0，迭代中止精度ε，最大迭代次数kmax，令迭代步k＝0；步骤2：对v0，计算目标函数值J和梯度g；设定初始搜索方向p0＝－g；令α0＝1／‖p‖2；步骤3：如果k＜kmax，迭代循环：1：用一维线搜索算法，得到搜索步长αk；2：速度模型修正vk＋1 ＝vk＋αkpk；3：对修正的vk＋1，计算新的目标函数值J和梯度g：4：若梯度或目标函数值满足终止条件，则算法终止，否则转步骤4；5：计算新的搜索方向pk；6：令k＝k＋1，返回步骤3继续.步骤4：输出结果v.在该算法中，需要计算目标函数的梯度g，对大规模模型的计算，我们基于下面的共轭方法来计算.2.3 梯度计算在BFGS方法中需要梯度的信息，这是一个M×（Nt×Nr）矩阵，这里Nr为接收点数.若用差分法求梯度需要求解M次的正问题，计算量太大.差分法计算梯度适合求解规模较小的问题（Zhang et al.，2013）.共轭法即求解原问题的共轭问题，是大规模问题中最常用的梯度计算方法，它只需要两次正问题的计算量就可以得到一个较为准确的梯度.下面给出这种方法的一般形式的推导.记则其中（xr，zr）为接收点位置坐标.对J（v）关于v微分，得：其中是u关于v的导数.因为u是关于v的非线性函数，不易直接求出，下面给出的方法是为了计算梯度时避免求解δu.为此，将正演问题写作由链式法则，得用任意的测试函数乘以上式且对时间和空间积分，得定义使得成立，所以：为消除δu，令即通过式（33）求出，则梯度可以由算出.若L为关于u的线性算子，则显然有对于波动方程（1），算子L关于u是线性的，且显而易见有下面考虑如何求的值.假设由于速度扰动δv引起的散射场的扰动为则v＋δv和u（v＋δv）同样满足声波方程：采用下列近似：将式（38）减去式（1），再将式（39）代入，且采用弱散射近似，可以得出方程：类似可得出与式（1）相同的初始条件.为求得算子，考虑到式（30），计算这里Δ是Laplace算子.考虑式（41）最后一个等号右端第1项，由Green公式可知：（42）式等号右端第2项中，已知且可令则第2项为0.等号右端第3项中，已知且可令则第3项为0.考虑式（41）最后一个等号右端第2项，由Green公式可知：注意到（43）式右端的边界项，在计算中希望边界项为0，因为δu的边界已知，所以根据情况选择使得边界项为0.要求上边界为零边界条件，其他三边使用吸收边界条件，显然零边界条件使得边界项为0，对于吸收边界条件，只要吸收得比较彻底，其实可以认为区域是无界的，这样就不需要考虑边界项的影响，可以看出好的吸收边界对于梯度的计算是很重要的，吸收不完全会对梯度计算的准确性带来影响.综合式（41）、（42）和（43），可以得到：由上述可知，共轭算子可描述为：关于“初始”时间的限制条件：边界条件与u相同，上边界为0，其他三边为吸收边界（8）—（10）.对于上述方程和定解条件，我们用时间逆推的方法进行求解，过程和计算量与求解正问题基本相同.最后给出梯度公式：3 数值计算本节用Marmousi复杂构造模型来进行反演数值计算.图1是Marmousi速度模型，常用来做偏移成像，用于验证各种成像方法的能力和效果（张文生，2009）.数值算例实验设计如下：计算的实际规模是地表宽度2952m，地下深度为1488m，记录时间为1.8s.空间步长Δx＝6m，Δz＝6m，由稳定性条件取Δt＝0.6ms.震源的最大频率为60Hz，中心频率为30Hz.离散后区域的空间点数为Nx＝493和Nz ＝249，时间取样点数为Nt＝3001.实验中共设80个炮点，每个炮点有40个接收点，双边接收.接收点之间的距离均为24m，炮检距离为24m.炮点和接收点置放在同一水平线上，深度为6m.图1 Marmousi速度结构模型Fig.1 Marmousi velocity structure model观测数据是由Marmousi模型正演得到的人工数据.图2是某炮的炮记录，图2a无边界吸收，图2b有边界吸收，比较可以看到，边界反射得到了明显的消除，这有助于提高梯度计算的精度.图3是用于反演的初始速度模型图，速度范围是1500m·s－1至4300m·s－1，中间每层线性递增得到的，每层速度相同；可以看到初始模型与真实模型的差异很大，已经完全没有原来真实模型的构造信息了. 频率多尺度方法就是根据问题对不同尺度的频率进行分解.为简单起见，我们采用截去某个频率以上的所有频率，只保留这个频率以下的频率进行反演，然后增加频率上界，依次进行反演.在本文计算中，选取的频率多尺度实验分为5个尺度：0～5Hz，0～10Hz，0～25Hz，0～35Hz，0～60Hz；每个尺度有逐级包含关系，这样可以更有效地恢复模型波数.为作为比较，在不采用本文频率多尺度方法的情况下也直接对原始数据反演，图4是迭代250步的结果，可以看到结果较差，再继续迭代也是如此.下面我们对5个频率尺度的波场进行逐级反演，每个尺度分别迭代50次.首先用0～5Hz频率尺度波场进行反演，图5a、5b和5c分别是迭代5步、20步和50步的反演结果，初始模型仍为图3.可以明显看到，反演结果随迭代次数增加逐步向好的方向收敛.为了节省篇幅，下面直接给出其他尺度迭代50次的反演结果.图6用是0～10Hz频率尺度的波场反演迭代50次的结果，初始模型是图5c；图7是用0～25Hz数据反演迭代50次的结果，初始模型是图6；图8是用0～35Hz数据反演迭代50次的结果，初始模型是图7；图9是用0～60Hz 数据反演迭代50次的结果，初始模型是图7.可以看到，反演结果从图6到图9逐步变好；图8和图9的结果比较接近，不过比较可以看出在细节上还是有些差别（比较下面的图10e和图10f更清楚）.为了进一步与真实模型作细致的比较，我们取CDP号为180处的位置，来与真实模型作比较；图10是5个不同频率尺度迭代50次的全波形反演结果（红线）与真实模型（蓝线）的比较.比较可知，在模型浅部，有非常好的收敛，在深部，反演的精度虽不如浅部，但仍能反映出良好的构造信息.原因有两方面，一是深部波场信息相对与浅层较弱，在计算中我们未对波场作任何处理；二是从数值计算的角度看，目标函数的未知量是一个向量，分量是每个网格点的速度值，目标函数梯度的分量对应每个网格点的梯度，从而也表示了不同深度的梯度，负梯度方向是未知量的修正方向，由于浅部接收信号强，浅部的变化量即梯度就会比较大，修正的时候着重于修正浅部.目标函数的二阶导数反映每个分量变化对梯度的影响，对深浅有平衡作用，但二阶导数是近似计算的，因此不可避免会导致深部反演结果的误差大于浅部.图11给出了0～50Hz频率尺度反演的目标函数值随迭代步的变化情况，可以看到前30次迭代，目标函数值下降很快，在80次之后，下降很小，表明继续迭代反演结果改善也不大，其他尺度的情况大致类似，考虑到计算量，每级反演设置最大迭代次数为50.此外，我们还对加噪声的情况也进行了数值计算，结果表明对20%的随机噪声也有校好的反演结果，为节省篇幅，不再详述.图2 某炮的炮记录（a）无边界吸收；（b）有边界吸收.Fig.2 A shot gather record with（a）no boundary absorption，and（b）boundary absorption 图3 初始速度模型Fig.3 Initial velocity model图4 不用频率多尺度方法，用0～60Hz的波场数据迭代250次的反演结果（初始模型为图3）Fig.4 The inversion result after 250iterations with 0～60Hz wavefield data.The frequency multiscale method is not adopted（The initial model is Fig.3）图5 0～5Hz的波场数据的全波形反演结果（初始模型为图3）（a）迭代5次；（b）迭代20次；（c）迭代50次.Fig.5 The full－waveform inversion results with 0～5Hz wavefield data（The initial model is Fig.3）（a）5iterations；（b）20iterations；（c）50iterations.图6 0～10Hz的波场数据迭代50次的全波形反演结果（初始模型为图5c）Fig.6 The full－waveform inversion result with 0～10Hz wavefield data after50iterations（The initial model is Fig.5c）The full－waveform inversion result with 0～25Hz wavefield data after50iterations（The initial model is Fig.6）图8 0～35Hz的波场数据迭代50次的全波形反演结果（初始模型为图7）Fig.8 The full－waveform inversion result with 0～35Hz wavefield data after50iterations（The initial model is Fig.7）最后指出，全波形反演是一个典型的大规模科学计算问题，计算量巨大，我们采用了MPI并行算法来实现.在计算中，对计算区域进行划分，分别赋予不同进程计算，这样，问题就分布存储在各个进程中，从而满足大规模计算对内存的要求.正问题是整个区域上关于时间递推的计算，每个离散点的计算都有它附近点进行参与，所以在进程管理区域边界与相邻的进程之间有信息交换.显然，边界点数与子区域点数比值越小，通信时间与计算时间的比值就越小计算效率就越高.通信过多会导致并行效率变差.由BFGS算法可知，搜索方向的计算只有向量参与，没有矩阵，向量点积只需在各进程中进行计算，再将所有进程所求的值进行相加即可，注意重叠区域不要进行多次运算.下面分析一下并行效率.我们知道，串行程序的执行时间近似等于程序指令执行花费的CPU时间，但并行程序相对复杂，其执行时间等于从并行程序开始执行，到所有进程执行完毕，墙上时钟走过的时间，也称之为墙上时间（Walltime）.对各个进程，墙上时间包括计算CPU时间、通信CPU时间、同步开销时间、同步导致的进程空闲时间.表1列出了4进程、8进程、16进程、32进程迭代反演10次的墙上计算时间的比较，由此可算出并行可扩展性.可以看到，平均的并行可扩展性大于0.9，是比较高的.计算在“科学与工程计算国家重点实验室”三号机群系统上完成，该机群有282个计算刀片，每个刀片包含两颗IntelX5550四核处理和24GB内存，单核双精度浮点峰值性能为10.68Gflops.用32个进程进行计算，总迭代250步，计算时间约47h.The full－waveform inversion result with 0～60Hz wavefield data after50iterations（The initial model is Fig.8）图11 0～50Hz频率多尺度反演目标函数随迭代步k的变化情况Fig.11 The variation of objective function vale of 0～50 Hz frequency multiscale inversion with the iterative number k表1 并行算法效率分析Table 1 Analysis of the efficiency of parallel algorithm 进程数墙上时间（s）并行可扩展性4 48094 8 24371 0.987 16 13418 0.89632 6656 0.903图10 在CDP180位置处，不同频率尺度波场迭代50次的全波形反演结果（红线）与真实模型（蓝线）的比较（a）初始模型；（b）0～5Hz；（c）0～10Hz；（d）0～25Hz；（e）0～35Hz；（f）0～60Hz.Fig.10 A comparison between the full－waveform inversion results after 50iterations on five different frequency scales（red line）and the exact model（blue line）at the position of CDP 180（a）Initial model；（b）0～5Hz；（c）0～10Hz；（d）0～25Hz；（e）0～35Hz；（f）0～60Hz.4 结论全波形反演是一个极小化模拟数据与已知记录之间残差的迭代过程.反演在时间域中进行，针对反演易陷入局部极小值的问题，我们发展了基于不同尺度频率数据的逐级反演方法，即用前一尺度的反演结果作为下一尺度反演的初值模型，有效地解决了初值离真解较远时，反演不收敛的问题.文中详细阐述和推导了理论公式和相应的算法，并基于Marmousi复杂构造模型进行了大规模反演计算，得到了较好的反演结果，说明了该方法的有效性和对初始模型具有较高的稳健性.波动方程全波形反演是一个典型的大规模科学计算问题，计算量巨大，文中用MPI并行方法来实现，提高了计算效率，这为进一步应用提供了基础.致谢本文的计算在科学与工程计算国家重点实验室三号机群系统上完成，在此表示感谢，此外还要感谢实验室很多老师的大力帮助和支持，在此不一一指出.最后，还要感谢审稿者提出的宝贵建议，使得本文能以目前的形式呈现.ReferencesBerenger J P.1994.A perfectly matched layer for the absorption of electromagnetic put.Phys.，114（2）：185－200.Bleistein N，Cohen J K，Hagin F G.1987.Two and one－half dimensional Born inversion with an arbitrary reference.Geophysics，52（1）：26－36.Broyden C G.1970.The convergence of a class of double rand minimization algorithm：2.The new algorithm.IMA J.Appl.Math.，6（3）：222－231. Bunks C，Saleck F，Zaleski S，et al.1995.Multiscale seismic waveform inversion.Geophysics，60（5）：1457－1473.Bursteddle C，Ghattas O.2009.Algorithmic strategies for full waveform inversion：1Dexperiments.Geophysics，74（6）：WCC37－WCC46.Choi Y，Min D J，Shin C.2008.Two－dimensional waveform inversion of multi－component data in acoustic－elastic coupled media.Geophysical Prospecting，56（6）：863－881.Clayton R，Engquist B.1977.Absorbing boundary conditions for acoustic and elastic wave equations.Bulletin of the Seismological Society of America，67（6）：1529－1540.Dahlen F A，Hung S H，Nolet G.2000.Frechet kernels for finite frequency traveltimes—Ⅰ.Theory.Geophys.J.Int.，141（1）：157－174.Davidon W C.1991.Variable metric method for minimization.SIAM J.Opt.，1（1）：1－17.Devancy A J.1984.Geophysical diffraction tomography.IEEE Transactions on Geoscience and Remote Sensing，22（1）：3－13.Epanomeritakis I，Akcelik V，Ghattas O，et al.2008.A Newton－CG method for large－scale three－dimensional elastic full waveform seismic inversion.Inverse Problems，24（3）：1－26.Fawcett J，Keller H B.1985.Three－dimension ray tracing and geophysical inversion in layered media.SIAM J.Appl.Math.，45（3）：491－501. Fletcher R，Powell M J D.1963.A rapidly convergent descent method for minimization.The Computer Journal，6（2）：163－168.Fletcher R.1970.A new approach to variable metric algorithms.The Computer Journal，13（3）：317－322.Gauthier O，Virieux J，Tarantola A.1986.Two－dimensional nonlinear inversion of seismic waveforms：Numerical results.Geophysics，51（7）：1387－1403.Gautier S，Nolet G，Virieux J.2008.Finite－frequency tomography in a crustal environment：application to the western part of the Gulf of Corinth.Geophys.Prospect.，56（4）：493－503.Goldfarb D.1970.A family of variable－metric methods derived by variational means.Mathematics of Computation，24（109）：23－26.Ha W，Chung W，Park E，et al.2012.2－D acoustic Laplace－domain waveform inversion of Marine field data.Geophys.J.Int.，190（1）：421－428.Hung S H，Shen Y，Chiao L Y.2004.Imaging seismic velocity structure beneath the Iceland hot spot：a finite frequency approach.J.Geophys.Res.，109（B8）：doi：10.1029／2003JB002889.Jin S，Madariaga R.1994.Nonlinear velocity inversion by a twostep Monte－Carlo method.Geophysics，59（4）：577－590.Mora P.1987. Nonlinear two－dimensional elastic inversion of multioffset seismic data.Geophysics，52（9）：1211－1228.Nocedal J，Yuan Y X.1993.Analysis of a self－scaling quasi－Newton method.Mathematical Programming，61（1－3）：19－37.Pratt R G，Worthington M H.1988.The application of diffraction tomography to cross－hole seismic data.Geophysics，53（10）：1284－1294.Pratt R G，Shin C，Hick G J.1998.Gauss－Newton and full Newton methods in frequency－space seismic waveform inversion.Geophys.J.Int.，133（2）：341－362.Pratt R G.1999.Seismic waveform inversion in the frequency domain PartⅠ：Theory and verification in a physical scale model.Geophysics，64（3）：888－901.Sambridge M，Drijkoningen G.1992.Genetic algorithms in seismic waveform inversion.Geophys.J.Int.，109（2）：323－342.Sen M K，Stoffa P L.1991.Nonlinear one－dimensional seismic waveform inversion using simulated annealing.Geophysics，56（10）：1624－1638. Sen M K，Stoffa P L.1995.Global Optimization Methods in Geophysical Inversion.The Netherlands：Elsevier Science B.V.Shanno D F.1970.Conditioning of quasi－Newton methods for function put.，24（111）：647－656.Sheen D H，Tuncy K，Baag C E，et al.2006.Time domain Gauss－Newton seismic waveform inversion in elastic media.Geophys.J.Int.，167（3）：1373－1384.Shin C，Cha Y H.2008.Waveform inversion in the Laplacedomain.Geophys.J，Int.，173（3）：922－931.Shin C，Ha W.2008.A comparison between the behavior of objective functions for waveform inversion in the frequency and Laplace domains.Geophysics，73（5）：VE119－VE133.Shin C，Cha Y H.2009.Waveform inversion in the Laplace－Fourier domain.Geophys.J.Int.，177（3）：1067－1079.Sirgue L，Pratt R G.2004.Efficient waveform inversion and imaging： A strategy for selecting temporal frequencies.Geophysics，69（1）：231－248.Song Z M，Williamson P R，Pratt R G.1995.Frequency－domain acoustic－wave modeling and inversion of crosshole data：Part 2—Inversion method，synthetic experiments and real－data results.Geophysics，60（3）：796－809.Stoffa P L，Sen M K.1991.Nonlinear multiparameter optimization using genetic algorithms：Inversion of plane－wave seismograms.Geophysics，56（11）：1974－1810.Tarantola A.1984.Inversion of seismic reflection data in the acoustic approximation.Geophysics，49（8）：1259－1266.Tarantola A.1986.A strategy for nonlinear elastic inversion of seismic reflection data.Geophysics，51（10）：1893－1903.Tong P，Zhao D P，Yang D H.2011.Tomography of the 1995 Kobe earthquake area：Comparison of ray and finite－frequency approaches.Geophys.J.Int.，187（1）：278－302.Yang D H，Liu E，Zhang Z J，et al.2002.Finite－difference modeling in two －dimensional anisotropic media using a fluxcorrected transport technique.Geophysical Journal International，148（2）：320－328.Yang D H，Wang S Q，Zhang Z J，et al.2003.N－times absorbing boundary conditions for compact finite difference modeling of acoustic and elastic wave propagation in the 2－D TI medium.Bulletin of the Seismological Society of America，93（6）：2389－2401.Yang W C.1993.Nonlinear Chaotic inversion of seismic traces，Part 1：Theory and numerical experiments.Chinese Journal of Geophysics，36（2）：222－232.Yang W C.1993.Nonlinear Chaotic inversion of seismic traces，Part 2：Lyapunov experiments and attractors.Chinese Journal of Geophysics，36（3）：376－386.Yuan Y X，Sun W Y.1997.Optimization Theory and Methods（in Chiese）.Beijing：Science Press.Zhang W S.2009.Wave Equation Imaging Methods and Computations（in Chinese）.Beijing：Science Press.Zhang W S，Luo J.2013.Full－waveform velocity inversion based on the acoustic wave equation.American Journal of Computational Mathematics，3：13－20.附中文参考文献袁亚湘，孙文瑜.1997.最优化理论与方法.北京：科学出版社. 张文生.2009.波动方程成像方法及其计算.北京：科学出版社.。

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

Open Journal of Natural Science 自然科学, 2020, 8(4), 258-263Published Online July 2020 in Hans. /journal/ojnshttps:///10.12677/ojns.2020.840342D Acoustic Wave Equation ForwardModeling in the Frequency DomainKun Han, Xiangchun Wang*School of Geophysics and Information Technology, China University of Geosciences (Beijing), BeijingReceived: Jun. 23rd, 2020; accepted: Jul. 6th, 2020; published: Jul. 13th, 2020AbstractForward modeling in frequency domain plays an important role in the numerical simulation of seismic waves. Compared with time domain forward modeling, frequency domain forward mod-eling has many advantages, such as suitable multi shot parallel operation, no time dispersion, flexible frequency band selection and small error. The coefficient matrix of different frequencies is relatively independent in the frequency domain forward modeling, which is suitable for the acce-leration of parallel computing and greatly improves the computing efficiency. In this paper, for the optimal 9-point difference scheme of frequency domain acoustic equation, the implicit expression and sparse matrix solution are studied, and the seismic wave field is simulated forward. The ac-curacy and validity of the method are verified by model calculation.KeywordsFrequency Domain, Forward Modeling, Acoustic Equation, Parallel Computing二维频率域声波方程正演模拟韩坤，王祥春*中国地质大学(北京)，地球物理与信息技术学院，北京收稿日期：2020年6月23日；录用日期：2020年7月6日；发布日期：2020年7月13日摘要频率域正演在地震波数值模拟中占有十分重要的地位。

相比于时间域正演，频率域正演具有适合多炮并*通讯作者。

韩坤，王祥春行运算、无时间频散、频带选取灵活、误差小等优点。

频率域正演时，不同频率的系数矩阵相对独立，适合并行计算加速处理，大大提高计算效率。

这里针对频率域声波方程的最优化9点差分格式，开展了隐式表达方式和稀疏矩阵的求解研究，正演模拟了地震波场。

模型计算验证了方法技术的准确性和有效性。

关键词频率域，正演模拟，声波方程，并行计算Copyright © 2020 by author(s) and Hans Publishers Inc.This work is licensed under the Creative Commons Attribution International License (CC BY 4.0)./licenses/by/4.0/1. 概述频率域波场正演相对于时间域数值模拟来说，有其自身的优势。

首先，在多炮数值模拟情况下，频率域相对于时间域效率更高，每个频率成分的阻抗矩阵只需要计算一次，加入并行计算后可以极大地提高计算效率；其次，就方程本身而言，频率域的衰减项更容易添加；另外，频率域模拟方法是计算区域所有网格点上的全部时间上的频率切片解，误差将分配到所有参与计算的每一个网格点上，因此不存在类似时间域模拟的累计误差，可以考虑长时间的地震波数值模拟。

Lysmer和Drake用有限元方法实现了地震波的数值模拟，初次提出了频率域正演模拟[1]。

Marfurt 指出频率域模拟地震波，能够非常有效地模拟出各向异性介质中波的传播及其衰减效应，而且适合大量炮点的同时模拟[2]。

Pratt和Worthington将频率域正演引入到波形反演中，应用常规二阶差分法(5点有限差分格式)对声波方程进行了离散化处理，结果表明当每个波长取至少13个网格点时，相速度误差就可以控制在1%以内，研究过程中详细推导了频率域声波方程和弹性波方程差分格式，构建了阻抗矩阵，指出频率域正演模拟适合多炮计算(每个频率成分共用一个阻抗矩阵)，便于模拟衰减效应，为频率域波场正演奠定了基础[3]。

由于频散严重和巨大的计算量等问题，频率域波场正演模拟没有得到广泛地应用，Jo和Shin等人在研究频率空间域声波波场正演模拟过程中，首次引入旋转坐标算子，将常规网格需要的5个网格点扩充到9个网格点，实现了频率域粘滞声波的正演模拟，频散得到了很好的压制，每个波长只需求4个网格就可以得到1%的精度要求，这使得频率域正演得到了极大地发展[4]。

Shin在Jo的基础上，由9点差分格式扩充到25点，相对提高了精度，这样在正演过程中单个波长内的网格点可以减少到2.5个，得到了理想的结果[5]。

进入21世纪后，频率域正演方法上有了更大的提高。

Min等人在Shin和Sohn研究的基础上做了进一步的研究，首先提出了一种新的25点有限差分加权最优化差分算子，通过Gauss-Newton法求出最优化加权系数，最后将频率空间域弹性波正演结果和方程的解析解进行了对比，从结果看出，利用新的差分算子计算的结果对于复杂介质来说结果比较理想，随后他又优化了25点算子的频散系数，通过最优化原理确定了新的加权系数，频散效应得到了更大限度的压制，并且减少了单位波长内所需的网格点数[6]。

高精度的时间域正演非常依赖于交错网格法，而频率域方法主要是基于混合网格法，Hustedt等人研究对比了频率域地震波模拟中交错网格法和混合网格法，指出对于二维模拟来说，交错网格法在内存需求和计算效率远不及混合网格法[7]。

Xu等从CPU计算时间、内存需求和硬盘读写速度三个方面研究了频率域声波模拟结果，并对比了LU分解和迭代法的优劣[8]。

韩坤，王祥春频域地震波数值模拟需要求解一个大型稀疏矩阵方程，该矩阵的带宽取决于正演算法(有限差分方法)，由于该矩阵具有高稀疏化、非对称性以及非正定性的特点，这些增加了求解的难度[9] [10]。

大型稀疏矩阵的求解方法分为直接法和迭代法。

直接法中应用较为普遍的是LU 分解方法，对于二维模型来说，LU 分解能够得到数值解，多炮情况下计算效率高，但是直接法对内存的需求较大。

迭代法对内存的需求小，但是并行计算效率不如直接法。

这里我们首先简要叙述了地震波频率域模拟方法的发展及研究现状，继而论述了有限差分数值模拟方法的基本原理，研究了频率域声波方程的最优化9点差分格式的隐式表达方式和稀疏矩阵的求解，最后通过正演算例验证了方法技术的准确性和有效性。

2. 频率域声波方程时间域声波方程为： ()()()()()2222222,,,,,,1,,,u x z t u x z t u x z t f x z t x z t v x z ∂∂∂+−=−∂∂∂式中的v 是介质的传播相速度，非均匀介质中为空间位置的函数，f 为震源函数，正演的时候通常选用雷克子波。

将上式变换到频率域，即对t 做傅里叶变换，(),,u x y t 变为(),,u x y ω，(),,f x y t 变为(),,f x y ω。

便有了下面的声波方程： ()()()()222222,,,,,,,,u x z u x y u x z f x z x z v ωωωωω∂∂+−=−∂∂ 3. 波动方程正演模拟有限差分数值模拟方法的原理是以离散代替连续，以差分代替微分，从而将波动方程差分化，这样就将一个微分问题转化为离散的数学问题[11] [12] [13]。

首先，将研究区域进行网格剖分；其次，采用有限差分算子将连续性方程在每个网格点处离散化；最后，求解差分离散后得到的差分方程组，得到每个网格点处解的近似值。

这里我们采用Jo 和Shin 等人提出的最优化9点差分格式进行计算。

根据最优化9点差分格式，可以构造出一个庞大的矩阵方程，该方程用矩阵运算简单表示为：()()(),A v P B ωωω=−。

式中，()B ω为模拟角频率为ω的雷克子波的傅里叶变换结果，矩阵()B ω只在震源处有值，其他处都为零。

其中(),A v ω是一个与模型参量、有限差分离散格式、频率和边界条件有关的矩阵，称为阻抗矩阵(Impedance Matrix)。

阻抗矩阵的大小取决于我们参与计算的网格大小，实际模拟过程中，阻抗矩阵往往是几万阶的矩阵，零元素占据大部分空间，所以我们称(),A v ω为大型稀疏矩阵。

求解此方程通常采用以下两种方法。

一种是直接法，此方法通常是对系数矩阵做LU 分解，进而用分解结果回代求解方程。

同频率的多炮正演共用一个分解结果，在二维情况下，计算效率较高。

系数矩阵虽然是一个稀疏矩阵，但是分解后的结果会有大量非零注入元，存储分解结果需要很大的存储空间，以致三维的存储变得不现实。

另一种方法为间接法，此方法通过迭代来求解方程的解，常用的方法有Jacobi 迭代、Gauss-Siedel 迭代、松弛迭代等；预条件对加速收敛的意义重大。

此方法避免了直接法的大矩阵存储问题，节约内存，但是计算效率偏低，失去了多炮共享分解结果的优势。

这里我们基于直接法进行求解。

4. 边界条件对于频率空间域地震波数值模拟来说，不设边界条件很容易产生明显的边界反射，从而对物理波场产生严重干扰，边界条件的构造非常重要，边界反射的压制效果直接影响着正演的精度。