二维频率域声波方程正演模拟
二维常密度声波波动方程
二维常密度声波波动方程
∂^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.波动方程成像方法及其计算.北京：科学出版社.。
FDMOD–声波方程有限差分正演模拟二维
FDMOD –声波方程有限差分正演模拟（二维）格式：fdmod <vfile >wfile nx= nz= tmax= xs= zs= [optional parameters] 必需的参数：<vfile 速度输入文件（包含速度值v[nx][nz]）>wfile 波场输出文件（包含每个时间步的波场值wave[nx][nz]）nx= x采样点个数(第二维)nz= z采样点个数(第一维)xs= 炮点x坐标dxs= 炮点x坐标间隔zs= 炮点z坐标dzs= 炮点z坐标间隔ns= 炮点个数tmax= 最大记录时间可选参数：nt=1+tmax/dt 时间采样点数(dt决定结果的稳定度)mt=1 波场输出时间切片的时间步长间隔dx=1.0 x采样间隔fx=0.0 x起始值dz=1.0 z采样间隔fz=0.0 z起始值fmax = vmin/(10.0*h) 震源子波的最高频率fpeak=0.5*fmax 雷克子波的峰值频率dfile= 密度输入文件（包含密度值d[nx][nz]）vsx= 垂直测线的x坐标hsz= 水平测线的z坐标rsx= 水平测线的起始检波器x坐标rlen= 水平测线长度rivl= 水平测线检波器采样间隔vsfile= 垂直测线的输出文件data[nz][nt]hsfile= 水平测线的输出文件data[nx][nt]ssfile= 震源点检波器的输出文件data[nt]verbose=0 =1 显示输出信息=2 更多输出信息abs=1,1,1,1 模型的顶,底,左,右使用吸收边界条件=0,1,1,1 顶部使用自由边界条件PML 参数:pml_max=1000.0 完全匹配边界条件参数pml_thick=0 完全匹配半层厚度(0 = 不使用PML边界条件)注意:程序使用普通显式二维差分方法。
二维弹性波频率域17点差分格式及正演模拟
二维弹性波频率域17点差分格式及正演模拟岳晓鹏;白超英;岳崇旺【摘要】为优化二维各向同性介质中弹性波频率域正演时阻抗矩阵的结构,减小正演所需内存,提高正演效率,在25点差分格式的基础上进行适当的简化,得到了二维弹性波频率域17点差分格式.该格式重新计算了弹性波中偏微分项和加速项的差分算子,减少了计算过程中的网格节点需求,构造了优化阻抗矩阵后的频率域正演矩阵方程;推导了纵波和横波相速度的频散公式,给出了不同泊松比条件下的频散曲线,得到了相速度误差控制范围±1％时每一横波波长内网格数需求.通过对比频散曲线和数值模拟时得到的波场快照及检波点处U、V分量,验证了17点差分格式与25点差分格式相比,具有稍严格的网格间距需求、相当的计算精度、略少的计算时间和更小的阻抗矩阵带宽等特点.%The 2D elastic wave frequency domain 17-point difference scheme is obtained based on the simplified 25-point difference scheme so as to optimize the structure of impedance matrix,reduce required memory and improve efficiency in 2D isotropic media forward.This scheme recalculates partial differential operators in the elastic wave equation,reduces the grid node in calculation process,and constructs the forward matrix equation with optimized impedance matrix.The P-wave and S-wave velocity dispersion formulae are derived,the dispersion curves of different Poisson's ratio are given,and the grid number requirement in each horizontal wave length is given with phase velocity error range being ± 1％.The frequency dispersion curve and numerical simulation of the wave field snapshots and seismic records show that the 17-point difference scheme has slightly strict grid spacing requirements,similaraccuracy,slightly less computation time and smaller impedance matrix bandwidth in comparison with 25-point difference scheme.【期刊名称】《物探与化探》【年(卷),期】2017(041)002【总页数】7页(P299-305)【关键词】弹性波方程;有限差分;频率域;数值模拟【作者】岳晓鹏;白超英;岳崇旺【作者单位】长安大学地质工程与测绘学院,陕西西安 710064;许昌学院数学与统计学院,河南许昌461000;长安大学地质工程与测绘学院,陕西西安 710064;长安大学地质工程与测绘学院,陕西西安 710064【正文语种】中文【中图分类】P631.4频率域正演模拟最早由Lysmer和Drake提出，他们利用该方法研究了多种介质中波的传播特征[1]；Shin等进一步发展了这种方法，将频率域有限元法应用于地震波形反演[2]；Jo等提出了最优化9点差分格式，极大地减小了数值频散效应[3]；Stekl和Pratt在Jo的基础上，将最优化9点差分格式引入到黏弹性介质的弹性波方程[4]；Min等提出了一种25点差分格式，通过计算最优化系数来减小数值频散，且减小了空间采样点数[5]；吴国忱教授分别将25点差分格式应用于VTI和TTI介质的正演模拟[6-7]；谷丙洛和梁光河等提出了一种21点差分格式对二维弹性波进行了数值模拟,该方法较25点差分格式能够保持精度相当且减少15%的计算内存和计算时间[8]。
声波波动方程正演模拟程序总结
声波波动方程正演模拟程序程序介绍：第一部分：加载震源，此处选用雷克子波当作震源。
编写震源程序后，我将输出的数据复制，然后我用excel做成了图片，以检验程序编写是否正确。
以下为雷克子波公式部分的程序：for(it=0;it<Nt;it++){t1=it*dt;t2=t1-t0;source[it]=(1.0-2.0*pow(PI*fm*t1,2.0))*exp(-pow(PI*fm*t1,2.0));fprintf(fp,"%.8f %.8f\n",t1,source[it]);}此处，为了成图完整，我用的是t2，而不是t1，也就是把雷克子波向右移动了一段距离，使主要部分都显示出来。
（频率采用的是30hz)从图中可以看出程序是正确的，符合理论上雷克子波的波形。
第二部分：主程序，编写声波正演模拟算子。
首先定义了各种变量，然后指定震源位置，选择权系数，给速度赋值，然后是差分算子的编写，这是主要部分，最后再进行时间转换，即把n-1时刻的值给n时刻，把n时刻的值给n+1时刻。
此处，我编写的是均匀介质声波方程规则网格的正演模拟程序，时间导数采用二阶中心差分、空间导数为2N阶差分精度，网格大小为200*200，总时间为400。
第三部分：这一部分就是记录文件。
首先记录Un文件，然后记录record文件。
模型构建与试算：1、我首先建立了一个均匀介质模型，首先利用不同时间，进行了数值模拟，得到波场快照如图所示：100ms 200ms 300ms此处，纵波速度为v=3000m/s。
模型大小为200×200，空间采样间隔为dx=dz=10m。
采用30Hz的雷克子波作为震源子波，时间采样间隔为1ms，图中可以看出，波场快照中的同相轴是圆形的，说明在均匀各向同性介质中，点源激发的波前面是一个圆，这与理论也是吻合的。
并且随着时间的增大，波前面的面积逐渐增大，说明地震波从震源中心向外传播。
声波方程正演模拟 (2)优秀课件
傅里叶变换法是利用空间的全部信息对波场函数 进行三角函数插值，能更加精确地模拟地震波的传播 规律，同时，利用快速傅里叶变换(FFT)进行计算，还 可以提高运算效率，其主要优点是精度高，占用内存 小，但缺点是计算速度较慢，对模型的适用性差，尤 其是不适应于速度横向变化剧烈的模型。
波动方程有限元法的做法是：将变分法用于单元分 析，得到单元矩阵，然后将单元矩阵总体求和得到总体 矩阵，最后求解总体矩阵得到波动方程的数值解；其主 要优点是理论上可适宜于任意地质体形态的模型，保证 复杂地层形态模拟的逼真性，达到很高的计算精度，但 有限元法的主要问题是占用内存和运算量均较大，不适 用于大规模模拟，因此该方法在地震波勘探中尚未得到 广泛地应用。
i 1, i, j 2 j 2
i 2, j2
x
图4-1 差分网格划分示意图
网格间隔长度
h
时间采，样步长
t
x ih
z jh
t kt
uk
uk i, j
i, j
表示(i,j)点k时刻的波场值
时间二阶、空间二阶差分格式推导如下：
将 uik,在j1 (i,j)点k时刻用Taylor展式展开：
u k 1 i, j
(4-8)
s(k为) 震源函数，
一般使用一个理论的雷克型子波代替，即：
e s(t) (2f / )2t2 cos 2ft
t 为时间
f 为中心频率，一般取为20-40HZ
为控制频带宽度的参数，一般取2-5
(i i0 ) * ( j j0 ) 确定震源位置
4.1 稳定性条件
对于特定的偏微分方程只有特定的几种有限差分 格式是无条件或有条件稳定的，（4-7）、（4-8）式即 是已被证明的有条件稳定格式，其稳定性条件分别为：
二维弹性波动方程组在频率域中的正演算法_郭华
第27卷 第3期 长春地质学院学报 V ol.27 N o.3 1997年7月 JOURNAL OF CHANGCHUN UNIVERSIT Y OF E ART H S CIENCES July 1997二维弹性波动方程组在频率域中的正演算法郭 华 马瑞杰(长春科技大学应用理学院,长春130026)摘要 运用有限差分法研究频率域中二维弹性波动方程组的正演算法及理论合成记录的制作。
频率域中二维弹性波动方程组的建立还未见报道,且此种算法要比时间域中正演算法计算量小得多,大大提高了计算速度,并且为弹性波动方程组的多参数反演提供了基本保证。
理论合成记录的制作是获取波场信息的重要途径,可以全面地反映地震波在地下介质中的分布与传播情况。
关键词 弹性波方程组 频率域 正演 合成记录 波场信息中图分类号 P 631.4第一作者简介 郭华 女 33岁 讲师 硕士 应用数学专业 已发表“频率域波动方程的双参数识别方法”等论文收稿日期 1996011 问题的提出波动方程描述波在介质中的运动和传播规律。
由于地下介质十分复杂,地震波型种类繁多。
近年来,弹性波动方程已开始直接应用于地震勘探的生产工作。
地震勘探的目的在于了解地球内部的结构、物性参数,由于它的不可直接测量的原因(或者测量的耗资太大)人们只好用能得到的某些信息去推断地下的结构和介质参数,因而这是数学中的反问题。
正演问题与反演问题是一个问题的二个方面,正演是基础,反演是目的。
2 数学模型的建立采用中间放炮的二维叠前地震记录进行研究。
由于地震波的传播满足弹性波运动的二阶偏微分方程组:(x i ) 2u i t 2= x i[ (x i ) u ]+∑3j =1 z [ (x i )( u i x j + u j x i )]+f i (i =1,2,3;j =1,2,3)其中:u =u 1i +u 2j +u 3k 为弹性体内质点的位移向量,x i ∈ ,f (x i ,t )是震源项, 、 是拉梅系数。
声波方程正演模拟共46页
t
x y z
vx 1 (P )
t x
vy 1 (P )
t y
vz 1 (P )
t z
11
二、波动方程类型及其局限性
能够描述且只能描述纵波的传播规律，包括 直达波、反射波、透射波、折射波等，但不能描 述转换波传播规律。
需要的已知条件包括： 1）震源函数 2）地层速度/密度 3）边界条件
对于二维速度-深度模型，地下介质中地震波 的传播规律可以近似地用声波方程描述：
2u t 2
v
2
(
2u x2
2u z 2
)
S (t )
(4-1)
v(x, z) 是介质在点（x , z）处的纵波速度，
u 为描述速度位或者压力的波场，
s(t) 为震源函数。
23
空间模型网格化（如图4-1所示）：
i 2, j2
12
2、弹性波方程：
2u t 2
v
p
2
(
2u x 2
2u z 2 ) S (t)
2
w
t 2
vs
2
(
2w x 2
2w z 2 )
vx
xx
xz
t
x
z
vz xz zz
t
x
z
xx ( 2u) vx vz
t
x
z
zz vx ( 2u) vz
t
x
z
xz u vz u vx
内容提纲ห้องสมุดไป่ตู้
一、地震勘探基本原理 二、波动方程类型及其 局限性 三、数值算法类型及其局限性 四、声波方程的有限差分法数值模拟
1
一、地震勘探基本原理
＊
毕设论文--粘声波正演模拟研究
本科毕业设计（论文）题目：粘声波正演模拟方法研究学生姓名：xxx学号：xxx专业班级：xxx指导教师：xxx2015年 6月20日粘声波正演模拟方法研究摘要地球上介质的黏滞性会引起大地的吸收效应,它会影响波场所有的频率成分,尤其对于高频的影响最大,导致地震分辨率降低。
黏滞吸收作用会影响地震波波形、频带、振幅等因素。
一个高效的粘声波正演模拟方法，可以考虑到由于实际介质造成的地震波的吸收衰减作用。
可以更加准确模拟地震波在非完全弹性实际地层中的传播，在这里，本文通过编程建立不同的粘声波方程数值模拟模型跟正常的声波方程数值模拟模型进行对比分析，从而了解粘声波正演模拟方法的优越性。
关键词：粘声波；正演模拟；有限差分；Study on the forward modeling of viscoelastic acousticwavesAbstractThe absorption effect is mainly caused by the viscosity of the earth media itself.The viscous stagnation can affect all the frequency components of the wave field.And the effect of the high frequency components is bigger,which leads to the decrease of seismic resolution.The absorption of the absorption has a great influence on the wave, frequency and amplitude of the seismic wave.. A highly effective viscoelastic forward modeling method can take into account the absorption and attenuation of seismic waves by real media.. Accurate simulation of the propagation of seismic waves in the actual strata of the imperfect elasticity. Here. In this paper, the program, establish different visco acoustic wave equation numerical simulation model with normal acoustic wave equation numerical simulation model for comparative analysis, to understand the visco acoustic forward modeling method of superiority.Keywords：Viscoelastic acoustic wave；Viscoelastic acoustic wave；Finite difference;目录第1章引言 (1)1.1 研究意义 (1)1.2 正演模拟方法 (1)1.3 国内外研究现状 (3)1.4 本论文研究内容 (3)第2章粘性介质基本理论 (5)2.1 粘性介质的基本特点 (5)2.2 粘性介质模型的构建 (6)2.2.1 开尔芬固体模型 (6)2.2.2 标准线性体 (7)2.3 品质因子 (7)2.3.1 定义 (8)2.3.2 与吸收系数 的关系 (8)2.3.3 与介质速度v的关系 (8)2.3.4 与震源频率的关系 (8)第3章二维各向同性介质粘声波方程数值模拟 (10)3.1 粘声波方程及其交错网格高阶差分格式 (10)3.1.1 粘声波方程的推导 (10)3.1.2 速度-应力方程的推导 (11)3.1.3 交错网格有限差分 (12)3.2 模拟震源 (16)3.2.1 震源的选择 (16)3.2.2 震源的类型 (17)3.3 边界条件 (18)3.3.1 衰减边界条件 (19)3.3.2 PML边界条件 (19)3.4 稳定性 (20)3.5 数值频散问题 (22)3.6 模型试算 (23)3.6.1地震波的传播规律研究 (23)3.6.2 粘声波的衰减规律研究 (27)第4章结论 (28)致谢 (30)参考文献 (31)第1章引言1.1 研究意义我国的石油天然气等资源类工业发展的非常快，随着油气资源的不断发现并且开采，现阶段所存在的油气资源越来越少，因此油气勘探工作变的越来越困难。
不同边界条件下的二维声波方程数值模拟
不同边界条件下的二维声波方程数值模拟姚铭;汪勇;杨晓柳;高刚;桂志先【摘要】地震波场数值模拟是研究波动现象的重要手段之一,对油气田的勘探和开发具有重要意义.数值模拟过程中,需要通过添加边界条件来尽可能消除由于截断所产生的边界反射.选取雷克子波作为震源项,分别建立均匀及层状地质模型,拟定合适的波场模拟参数,实现了不同边界条件下的二维声波方程数值模拟.利用数值模拟得到的波场快照和地震记录直观地对比分析不同边界条件对边界反射的消除效果,认为透明边界条件(TBC吸收边界条件)和Clayton-Engquist边界条件(CE吸收边界条件)都能够较好地消除边界反射.最后提出了一种组合边界条件的方法.%Numerical simulation of seismic wave field is one of the important means to study the fluctuation phe-nomenon , which is very important for the exploration and development of oil and gas field .In the numerical simula-tion process, it is necessary to eliminate the boundary reflection due to the truncation by adding the boundary condi -tion as much as possible .The Ricker is selected as the source , and the uniform and stratiform geologic model is es-tablished respectively .The appropriate wave field simulation parameters are developed to realize the numerical sim-ulation of the two-dimensional acoustic equation under different boundary conditions .Wavelet snapshots and seismic records obtained by numerical simulation are used to visually analyze the elimination effect of boundary conditions on boundary reflections .It is considered that the boundary condition is better to eliminate the boundary condition( hereinafter referred to as TBC absorption boundary condition ) andClayton-Engquist boundary condition ( hereinaf-ter referred to as CE absorption boundary condition ) .Finally, a method of combining boundary conditions is pres-ented.【期刊名称】《科学技术与工程》【年(卷),期】2017(017)032【总页数】6页(P11-16)【关键词】边界条件;声波方程;数值模拟;有限差分【作者】姚铭;汪勇;杨晓柳;高刚;桂志先【作者单位】长江大学油气资源与勘探技术教育部重点实验室,地球物理与石油资源学院,武汉 430100;长江大学油气资源与勘探技术教育部重点实验室,地球物理与石油资源学院,武汉 430100;长江大学油气资源与勘探技术教育部重点实验室,地球物理与石油资源学院,武汉 430100;长江大学油气资源与勘探技术教育部重点实验室,地球物理与石油资源学院,武汉 430100;长江大学油气资源与勘探技术教育部重点实验室,地球物理与石油资源学院,武汉 430100【正文语种】中文【中图分类】P631随着计算机技术的快速发展，数值模拟方法不仅被广泛应用于勘探地震学的研究，在地震资料采集、处理和解释的各个环节中也都离不开它，它对于石油等能源的勘探和开发具有重大的意义。
《物探化探计算技术》 2019年1~6期总要目
物探化探计算技术
Ｖｏｌ．４１ Ｎｏ．６
２０１９ 年 １１ 月
ＣＯＭＰＵＴＩＮＧ ＴＥＣＨＮＩＱＵＥＳＦＯＲ ＧＥＯＰＨＹＳＩＣＡＬ ＡＮＤ ＧＥＯＣＨＥＭＩＣＡＬＥＸＰＬＯＲＡＴＩＯＮ
Ｎｏｖ．２０１９
Байду номын сангаас
《物探化探计算技术》
２０１９年１～６期总要目
第１期
一种适用于欠压实成因的地层孔隙流体压力预测技术 熊晓军，黄 劲，陈 容，ＬＩＡＯ Ｙｉｄｕｏ，袁 野（ １ ） 基于 Ｌ１范数的波阻抗反演 刘百红，李建华，郑四连（ ６ ） 针对岩性圈闭的正演照明采集参数优化设计 罗 伟，阎 贫，万琼华，陈兆明，黄 鑫（１２ ） 基于 ＬＳＱＲ 算法的二维声波方程频率域正演模拟与数值实现 张鑫磊，陈建宇（１８ ） 渤海 Ｑ 油田开发中后期高精度变速成图方法研究 陈华靖，张鹏志，任百聪，贾海良（２７ ） 岩石物理建模技术在致密砂岩储层预测中的应用———以鄂尔多斯盆地北部 Ｈ 区块为例 王震宇，刘俊州（３４ ） 强反射能量屏蔽补偿在渤海 Ｃ 构造中的应用 周 星，徐德奎，何 玉（４１ ） 煤层瓦斯隧道中开展 ＴＳＰ 探测工作的研究与探索 智 刚（４７ ） 基于解析雅克比矩阵的 Ｅ－Ｅｘ广域电磁法一维并行约束反演 索光运，李帝铨，胡艳芳（５５ ） 垂直接触带上电阻率异常的空间分布特征 李祖强，熊 彬，兰怀慷，罗天涯，梁 卓（６２ ） 不同噪声水平高密度电法的分辨率和勘探深度研究 邢润林，陈儒军，刘海飞，王小杰，陈兴生（６８ ） 三维地质体对 ＡＭＴ 二维反演的影响研究 张逗逗，汤井田，肖 晓，王 宁（８０ ） ＸＧＹ－２０２０Ａ 型原子荧光光谱仪的研制 李 可，张 勤（９１ ） 场源不同近似处理在直流电法边界元正演中的应用研究 梁 卓，熊 彬，徐志峰，兰怀慷，
声波方程数值模拟实验报告
声波方程数值模拟实验报告学号：201107010230姓名：包灵指导老师：程冰洁老师完成时间：2014年11月26日星期三实验要求：1、应用声波方程作为正演模拟的波动方程；2、将所提供震源函数离散后绘图；3、给定两个二维速度-深度模型（一个小模型；一个大模型），绘出图形来；4、对于小模型，整个区域的速度值可设为常数，即只有一种介质，将震源点放在模型中间，分别记录两个时刻的波前快照（即该时刻区域内所有网格点的波场值）。
第一时刻为地震波还未传播到边界上的某时刻，第二时刻为地震波已经传播到边界上的某时刻，体会其人工边界反射；5、对于大模型，定义为水平层状速度模型（至少两层）；做两个实验，一是将震源点放在区域表层任一点，记录下某些时刻的波前快照，体会地震波在两种介质的分界面上传播规律；二是合成一个地震记录，即记录下与震源同一深度点的各点所有时刻的波场值，并指出记录上的同向轴分别对应哪些波。
实验目的：1.通过本次作业，加深对波动方程的理解，明白波动方程所代表的物理意义。
2. 通过模拟地震波在介质中的传播，理解实际勘探中地震波在地层中的传播规律。
3. 通过模拟水平层状速度模型，体会地震波在两种介质分界面的传播规律，并能够从地震记录中识别出反射波，透射波，多次波，折射波和绕射波。
4. 通过模拟人工合成的地震记录，体会地震勘探基本原理和方法，验证地震波传播能量波形变化趋势。
需要的已知条件包括：1）震源函数2）地层速度（波速）3）边界条件2．弹性波方程：⎪⎪⎩⎪⎪⎨⎧∂∂+∂∂=∂∂+∂∂+∂∂=∂∂)()()(22222222222222z w x w v t w t S zu x uv t u s p声波方程的有限差分法数值模拟对于二维速度-深度模型，地下介质中地震波的传播规律可以近似地用声波方程描述：)()(2222222t S zu x u v t u +∂∂+∂∂=∂∂ （4-1） (,)v x z 是介质在点（x , z ）处的纵波速度，u 为描述速度位或者压力的波场，)(t s 为震源函数。
- 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. 边界条件对于频率空间域地震波数值模拟来说，不设边界条件很容易产生明显的边界反射，从而对物理波场产生严重干扰，边界条件的构造非常重要，边界反射的压制效果直接影响着正演的精度。