第 1 章 引言
1.1 问题背景
石油是工业的血液，其重要性不言而喻。随着我国经济的进一步发展，对石油、天然气等能源的需求也越来越迫切，大力提高国内相关行业的油气资源勘探水平成为当务之急。作为地球演化的产物，石油的生成、贮藏以及在地下的运移与地质结构密切相关。更好地认识地球内部构造可以大致知道石油的分布和储存情况，从而更好地适应实际勘探需要，更好地指导石油工业生产。目前，地质构造的研究主要采用地质勘探方法，如地面地质调查、各种物探方法(如重力勘探，磁力勘探，电法勘探、地震勘探等)，由于地震勘探准确、方便、经济、有效，所以地震勘探方法在地质勘探中应用最为广泛。
地震勘探方法采用人工爆破法、落重法等方式激发地震波，通过研究布置在地表的台站所接收到的地震信号，反推地下岩层构造，进而找出油气、煤矿或其他矿产的分布位置。从数学上讲，这是一个波动方程反演的过程。目前地学中常用的反演方法(如层析成像、地震偏移等)都是建立在正演的基础之上的，所以研究一种快速、高效和高保真的正演技术具有重要的意义。
1.2 地球介质与弹性波方程
弹性波方程是弹性介质中波的传播所满足的偏微分方程，也是当前地震勘探领域研究最多的方程之一。在单相弹性介质中，弹性波方程由以下方程构成：
由牛顿第二定律，可知当弹性介质在应力作用下处于运动状态时作用在物体上的面力、体力的合力与惯性力的合力为零，从而得到弹性介质运动平衡方程(位移——应力关系)：
完全的各向异性介质(三斜对称各向异性介质)中的弹性波方程有 21 个独立的弹性参数，非常复杂，很难研究弹性波在这种完全各向异性介质中的传播规律。同时，由于地层介质具有各式各样的对称性，利用这些对称性，可以大大简化弹性波方程，有利于更好地认识弹性波在地层中的传播规律。各向同性介质和横向各向同性介质是地学中研究做多的两种介质。下面分别介绍在这两种介质中传播的弹性波方程。
1.2.1 各向同性介质中的弹性波方程
在各向同性介质中，弹性波沿着各个方向的传播速度相同，四阶弹性张量中只有两个独立参数，弹性张量矩阵可以写为
第 2 章 三维 WNAD 算法及其频散分析
与传统高阶有限差分不同的是，NAD 类算法不仅使用位移，同时使用位移的梯度共同重构波场。梯度可以反映函数的变化趋势，梯度信息的使用使得 NAD类算法即使每个最小主波长内使用较少的采样点，仍然可以得到低数值频散、高分辨度的地震波场。本章将给出三维 WNAD 算法的计算格式并分析其稳定性，相速度和群速度的频散分析表明 WNAD 算法较传统方法能更好地压制波场模拟中的数值频散。
2.1 加权的近似解析离散化（WNAD）方法
数值频散是差分方法数值离散控制方程而产生的固有现象。传统的差分方法在进行大规模地震波场模拟时很容易陷入两难境地：一方面，为了准确刻画一个波形，得到低数值频散的波场记录，需要在每个最小主波长内使用较多的网格点；另一方面，在计算内存和计算速度有限的前提下，为了进行大规模地震波场模拟，希望能够使用尽可能大的空间步长以减小计算所需的存储空间，同时加快计算速度，而采用大的空间步长意味着每个最小主波长内需要使用尽可能少的网格点数。这是一对矛盾。认真分析传统高阶差分格式，不难发现传统的一般高阶差分格式(如图(2.1a)所示)仅仅使用插值节点处的位移来重构插值节点处的高阶空间偏导数，波场信息中只有位移信息在下一步计算中得以保留，这使得一般高阶差分方法每个最小主波长内需要使用较多的网格点以消除数值频散。鉴于此，Yang 和他的合作者在充分研究了波场重构理论后，认为梯度反应了函数变化趋势，因此梯度信息非常重要。在重构波场时应充分考虑梯度的作用，在差分格式的构造中利用梯度、位移共同重构下一时间层的波场，而不是像传统差分方法那样对梯度信息弃之不用。求解波动方程的近似解析离散化类(NAD-type)方法就是在这种思想的指导下被提出并发展起来的。由于 NAD 类算法(如图 2.1b)使用位移和梯度共同重构插值节点的高阶偏导数，这使得这类方法可以使用较少的网格点得到高分辨率的、低数值频散的波场记录，在粗网格(大空间步长)条件下可以有效压制弹性波场中的数值频散，具有重要的理论和实际意义。
加权近似解析离散化(WNAD)方法是近似解析离散化类方法的一种，目前 WNAD 算法已经被应用于二维弹性波场模拟。本章首先将二维 WNAD 算法推广到三维情形，得到三维 WNAD 算法，然后进行相关的理论分析(如稳定性分析、频散分析等)。
2.2 WNAD算法的频散关系
物理上均匀各向同性弹性介质中传播的弹性波和声波的速度是一个与频率无关的常数。但是，因为数值求解声波、弹性波方程时，波传播方程和差分方程之间不可避免地存在误差，这些误差使得数值波速与频率相关，即产生所谓的“数值频散”现象。这些数值频散是客观上不存在的，它们严重干扰了地学工作者识别偏移等结果中的有效信号，甚至会导致虚假的地层构造解释，对实际生产产生错误的指导，所以需要选择一种能有效压制数值频散的波动方程数值求解方法。在使用数值方法进行弹性波正反演计算之前，需要知道这种数值算法得到的数值结果是否真实可靠。频散关系分析是判断数值方法优劣的重要方法，下文将逐次研究一维、二维和三维的 WNAD 算法的频散关系。
一维波动方程的频散关系是研究高维频散关系的基础，所以本文首先给出并比较求解一维波动方程的中心差分(FD)方法、四阶 Lax-Wendroff 修正 (LWC) 格式和 WNAD 算法的数值频散关系。
2.3.1 一维波动方程的频散分析
考虑如下的一维波动方程：
第 3 章 并行的加权近似解析离散化方法...................42
3.1 并行计算简介 ............... 42
3.1.1 并行计算平台.............. 42
3.1.2 并行程序设计简介............. 43
第 4 章 三维横向各向同性介质中的弹性波场模拟..............60
4.1 三维均匀各向同性介质模型 ................. 61
4.2 三维均匀 VTI 介质 ............... 66
第 5 章 WNAD 算法中的水平自由地表条件处理 ....................84
5.1 弹性波方程和自由地表条件 ............................. 84
5.1.1 三维弹性波方程及自由地表边界条件.................. 84
第6 章 基于 WNAD 的非一致网格算法
通过上一章的数值算例，不难看出WNAD算法在处理自由地表附近需要使用较高的空间采样率才能有效压制数值频散，但在计算区域内部使用相对较低的采样率即可有效压制波场中的数值频散。另外，一般来说，在介质变化剧烈区域(如含薄互层、低速区、含盐丘类高速区等复杂地质构造)，差分算法也需要采用较高的空间采样率才能有效压制数值频散，但是在介质变化缓慢的区域则不需要使用高的空间采样率。据此，本章研究并提出基于 WNAD 算法的非一致网格算法，在不同的计算区域采用不同的采样率，在能够有效压制数值频散的基础上，提高算法的计算效率，降低算法的存储需求。
6.1 非一致网格算法
有限差分格式在低速区或者特殊构造区域(如自由地表附近、内部复杂物性界面等)需要采用细网格以压制数值频散、提高计算结果的分辨度。研究弹性波在含低速区或者特殊构造区内的传播时，如果仅使用粗网格进行空间划分，则会产生数值频散，并降低整个波场的分辨度；而使用细网格划分，在介质变化平缓地区、高速区采样率会高于压制数值频散所需要的最低采样率，内存和计算存在一定的浪费。因此，只有根据地质模型特点，在不同的区域应当采用不同的采样率，这样才能兼顾计算效率和压制数值频散，所以需要研究并使用非一致网格算法进行弹性波场模拟。
非一致网格算法的基本思想是：根据介质速度大小(高速区采用粗网格，低速区采用细网格)或实际计算的需要(边界介质突变区域采用细网格，其他区域采用粗网格)等，在不同的区域选择不同的离散网格步长。1989 年，Moczo首次提出了非一致网格算法的思想，并进行了二维SH地震波场模拟，Jastram应用非一致网格算法模拟了弹性波场。国内学者如孙卫涛、李胜军、朱生旺、黄超、宋国杰等也相继将非一致网格剖分方法应用到了地震波场的数值模拟中。
本节将给出基于加权近似解析离散化(WNAD)方法的非一致网格算法，以期在保持能够有效压制数值频散的前提条件下，实现计算效率和存储效率的进一步提高。非一致网格算法实现的关键在于解决好网格变化过渡区域的联结问题。保证在过渡区域不因网格尺度变化而产生明显的人工绕射或者数值频散等人工噪声。
第7 章 结论、认识与展望
7.1 结论和认识
本文将求解声波、弹性波方程的 WNAD 算法推广到三维情况，获得了三维WNAD 算法，并发展了并行化的三维 WNAD 算法。在此基础上，论文研究了三维WNAD 算法的稳定性，分析了一维、二维、三维 WNAD 算法的频散关系，给出了并行WNAD算法的并行加速比和并行效率。结合应力镜像法，给出了WNAD算法中的自由地表条件及其离散化处理方法，并针对模拟中出现的新问题提出了基于WNAD算法的非一致网格算法，进一步提高了地震波场的模拟效率。本文的主要成果、结论及认识如下：
1. 发展了三维加权近似解析离散化(WNAD)方法，并分析了三维 WNAD 算法的稳定性。
2. 研究了WNAD算法的频散关系，从理论分析的角度阐明了WNAD算法较传统的有限差分格式能够更有效地压制波场中的数值频散，且WNAD算法的数值频散各向异性较弱，能够得到更加准确的波形。以三维情形为例，当每个最小主波长内的空间采样点数为2时，对于沿着坐标轴方向传播的地震波，WNAD算法得到相速度和群速度误差分别为8%和40%，而LWC方法的相速度和群速度的误差达到25%和100%，二阶FD差分格式的相速度和群速度误差达到36%和100%；在体对角线上，WNAD算法的相速度、群速度误差分别为3%和10%，而四阶LWC方法的误差达到了10%和25%，二阶FD算法的误差达到了16%和40%。数值频散分析表明，在同样的网格条件下，WNAD算法可以得到更加准确的高频信号；同时也表明，在使用大步长进行波场模拟时，WNAD算法能够模拟获得更为准确的地震波相位和振幅，可以有效压制波场中的数值频散。
3. 提出了并行化的WNAD算法，分析了并行WNAD算法的并行加速比和并行效率。和并行的SG方法相比，并行WNAD算法和并行SG方法的并行加速比均随着参与计算的进程数增多而增大，但是并行WNAD算法的加速比与进程数的比值大约为0.5，而并行SG方法的并行加速比和进程数的比值大约为0.1，说明并行WNAD算法的并行加速比和并行效率均优于SG方法。这是因为：(1)差分格式的并行效率和并行加速比与差分算子长度成反比，WNAD算法的算子长度仅为3，算子的局部性好；(2)WNAD算法可以使用大步长进行地震波场模拟，所以WNAD算法的网格点数比SG方法的网格点数少，这进一步降低了并行WNAD算法的通信量。事实上，在使用同样多的进程数进行计算时，要达到同样的消除数值频散效果，并行WNAD算法的计算速度是并行SG方法的6.38倍，通信量只需并行SG方法通信量的14.8%，存储量也只需并行SG方法存储量的14.8%。这说明并行WNAD算法可以使用较少的进程数得到更高的并行加速比，对并行计算机的利用率更高。
参考文献（略）