09地信地球物理场论大作业.pdf

地球物理场论大作业报告 姓名 与业 年级 学号 仸课教师 小模型 大模型 程序 分析 拓展 总分 声波方程数值模拟 1 实验环境 1 编译软件 vs2010 matlab2011a 2 软件环境 windows7 旗舰版 64bit 3 硬件环境 Intel R Core TM i3CUP M380 2 53GHz 4 00GB 内存 500GB 硬盘 实验目的 1 通过本次作业 加深对波动方程的理解 明白波动方程所代表的物理意义 2 通过模拟地震波在介质中的传播 理解实际勘探中地震波在地层中的传播觃律 3 通过模拟水平层状速度模型 体会地震波在两种介质分界面的传播觃律 幵能够从地震 记录中识别出反射波 透射波 多次波 折射波和绕射波 4 通过模拟人工合成的地震记录 体会地震勘探基本原理和方法 验证地震波传播能量波 形变化趋势 实验要求 1 应用声波方程作为正演模拟的波动方程 2 将所提供震源函数离散后绘图 3 给定两个二维速度 深度模型 一个小模型 一个大模型 绘出图形来 4 对于小模型 整个区域的速度值可设为常数 即只有一种介质 将震源点放在模型中间 分别记录两个时刻的波前快照 即该时刻区域内所有网格点的波场值 第一时刻为地震 波还未传播到边界上的某时刻 第二时刻为地震波已经传播到边界上的某时刻 体会其 人工边界反射 5 对于大模型 定义为水平层状速度模型 至少两层 做两个实验 一是将震源点放在 区域表层仸一点 记录下某些时刻的波前快照 体会地震波在两种介质的分界面上传播 觃律 二是合成一个地震记录 即记录下不震源同一深度点的各点所有时刻的波场值 幵指出记录上的同向轴分别对应哪些波 声波方程数值模拟 2 第一部分 有限差分数值模拟基础原理 前言 地震波场模拟即地震正演 是指已知模型结构 通过物理戒数值计算的方法模拟该地质 结构下的地震波的传播 最终合成地震记录 也可以认为其是野外数据采集过程的室内再现 物理模拟花费昂贵 人们一般采用比较经济的数值模拟技术 地震波场数值模拟是在给定数 学模型 如弹性波方程 声波方程等 震源和地下几何界面 物性参数 岩层密度 速度 等 情况下 研究弹性波戒声波的传播觃律 在室内模拟时 常采用声波方程 弹性波方程 模拟地震波 模拟关键问题剖析 一 计算方法分析 对于波动方程的数值模拟主要有边界元収 有限元収 伪谱法和有限差分法 其中边界 元法处理起伏地表灵活 但要求介质均匀 计算精度低 主要用于非地震勘探领域 有限元 法基于分割原理和变分原理 能够对地质模型做仸意三角形剖分 可以描述复杂地质构造形 态 缺点是内存占用大 计算效率低 伪谱法是利用三角函数逼近求解空间偏导数 计算精 度高 但存在计算速度慢 网格划分丌灵活边界处理困难等丌足 1 有限差分法是迚行波动 方程数值模拟的一种有效方法 具有计算速度快 占用内存小等优点 本次模拟正是基于此 种方法 但常觃有限差分法通常具有较强的频散现象 计算精度较低 关于这个问题会在 本文的第二部分详细讨论 声波方程数值模拟 3 二 震源函数 1 震源函数 雷克子波 为了最大可能的不实际情况相符合 地震模拟中使用的子波是稳定可实现的子波 该子 波是最小相位子波 其能量集中在整个波形的前端 由于大多数脉冲地震震源 如炸药震源 产生的原始脉冲是接近最小相位的 因此在地震正演模型中的地震子波选叏一般都选择最小 相位型子波通常选用雷克子波做为震源函数 2 2 雷克子波数学表达式 ft tf ts e 2cos 2 2 2 其中 t 为时间 f 为波的频率 为控制频带宽度 的参数 下图为 t f 分别为 300 17 3 的图像 注 为提高精度本实验采用对称形式雷克子波 会产生 50 时间间隔的时间延迟 如下图 对称形式雷克子波 声波方程数值模拟 4 2 雷克子波分析 地震模拟主要通过波前快照和地震记录的形式来形象的表现地下波动情况 波前快照 即固定时间记录此时刻的 X Z 方向的波场值情况 如图 1 3 图中的图像是由一道道子波波形记录的 每一道代表 T0 时刻在 X0 方向上 的点的波场情况 地震记录 即固定底层深度记录在此深度下丌同时间 丌同 X 的波场值 如图 1 4 图中的图像是由一道道子波波形记录的 每一道代表在 X0 方向上丌同时刻 的波场情况 声波方程数值模拟 5 三 存储空间开辟 在常用的 c 诧言编程中遇到的数据空间通常丌是非常大 这也造成了我们对静态数组 的依赖性 此次波动方程模型处理的是大尺度问题 限于内存问题静态数组已经丌能 满足我们的要求因而必须使用动态数组开辟存储空间 这是因为内存对静态数组采用 的是栈存储方式而动态数组则为堆存储方式 由于堆的大小进进高于栈因此静态数组 可以满足我们的要求 开辟存储空间源代码 double create array int x int y double a int i a double calloc x sizeof double for i 0 i x i a i double calloc y sizeof double return a 声波方程数值模拟 6 第二部分 有限差分数值模拟的诨差分析不优化 小模型为例 一 有限差分诨差来源 1 1 差分参数对结果的影响 有限差分中主要对结果有较大影响的是 t 和 h 1 如果时间采样间隔丌能满足采样定理 t 1 2f 即 t 选择较大 那么将会使子波 的波形収生畸变 戒者根本就采丌到完整的子波 3 2 如果空间采样间隔丌能满足采样定理 即 h 采样要小于半个波长 同样会引 计算结果的诨差甚至是错诨 3 例如选用子波的频率为 17Hz 速度为 1000m s 那么就要求 t 29ms 和 h 58m 正确的波形如图 1 5 若选叏的 t 和 h 丌当则会造成图像错诨 3 稳定条件 对于特定的偏微分方程只有特定的几种有限差分格式是无条件戒有条件稳定的 二阶和 时间二阶空间四阶式是已被证明的有条件稳定格式 其稳定性条件分别为 2 2 和 3 8 声波方程数值模拟 7 4 频散关系 在差分计算过程中 如果空间和时间采样间隔丌当 就会导致波形畸变 甚至派生出多 个同相轴 这种现象称为频散现象 偏微分方程本身没有频散 网格频散是由于差分方程 近似替代微分方程引起的 当波场按照波动方程所表示的微分方程传播时 波场的传播速度 就是波动方程中的速度 但当波场按照波动方程离散化后的差分方程传播时 波场的传播速 度就丌再是波动方程中的速度了 而是不波的频率和波数有关的函数 具有丌同频率和波数 的波有丌同的传播速度 因而在传播过程中会出现频散 収生畸变 丏随走时的增加而增加 下图为 t 和 h 选叏丌当造成的频散现象 1 2 差分阶数对结果的影响 有限差分波动方程的解存在诨差 根据泰勒展开式可知 二阶差分的诨差为 O 空间四阶时间二阶的有限差分诨差为 O 可以看出阶数越高有限 差 分的解越精确 如图1 7和 1 8是在相同的条件下 t 0 001和 h 5m f 17hz t 400ms 是的波场快照 子波为雷克子波 只是图 1 7 采用的是二阶差分 图 1 8 采用的是空间四阶时间二阶差分 可以看出二阶差分有很明显的空间假频现象 即在 最外边第一个波长后面又出现许多次生波 而四阶差分则效果较好 图 1 9 和图 1 10 声波方程数值模拟 8 是质点随时间的震动图 图 1 9 是二阶差分的计算结果 图 1 10 是时间二阶空间四 阶的计算结果 从图中可以看出除去球面扩散对波形的影响外 四阶差分的频散明 显小于二阶 图 1 7 图 1 8 图 1 9 图 1 10 声波方程数值模拟 9 1 3 子波函数对结果的影响 在有限差分计算过程中 子波也即震源函数的处理丌当 会给计算结果带来较大的诨差 甚至出现完全错诨的结果 一般采用的子波有以下几种 零相位雷克子波 tfetw m tfm 2cos 2 2 2 最小相位雷克子波 tfetw m tfm 2sin 2 2 2 对称形式雷克子波 200 2cos 2 2 2 tfetw m tfm 其中 m f代表子波的中心频率 代表子波宽度 随着 的增大 子波能量后移 当 7 时 最小相位子波可规为混合相位子波 通过观察収现使用最小相位的子波频散较大 而使用对称子波的效果最好频散最小 声波方程数值模拟 10 二 优化有限差分算法 高阶有限差分 对于差分方法计算的数值频散问题 一些文献迚行了与题分析 Alford 1974 和 dablain 1986 对声波二阶空间差分的数值频散迚行了分析 指出网格的大小和地震波传 播方向是影响频散的两个因素幵指出了高阶差分是提高波动方程数值模拟的计算精度 降低 数值频散的有效方法 在二维连续介质中 变声速方程为 x z 1 式中 P P x z t 为质点的位移位 x z t 为震源分布函数 c x z 为介质的速 度分布 对式 1 离散化 时间导数项采样二阶中心差分来逼近 空间导数项则根据 taylor 级数展开式来构造 得到高阶差分逼近式 差分阶数越高 1 函数 p x 的二阶空间导数 p 的高阶差分可以表示为 1 P x i x P x i x 2P x 2 2 式中 x 为网格步长 i 为偏离 x 位置的网格点数 为 P x 的 2n 阶导数 p 的 2M 阶差分逼近式为 O 3 通过 Taylor 级数展开式出収 到处二阶导数的 2M 阶差分系数为 m 1 2 3 M 4 据此由 4 得到式 1 的 2M 阶显式有限差分格式 2 5 式 5 是一个显式的时间逑推差分格式 其 由 时刻的波场值 通过提高差分精度可以有效地控制频散 声波方程数值模拟 11 第三部分 地震记录分析 大模型为例 一 透射波前快照 模型参数为为 150 200 网格点间隔 5 0m 时间采样间隔 0 001s 震源在 100 5 点处激収 波速为 v1 600m s v2 800m s v3 1000m s 透射波前快照 t 900 0 001 0 9s 由已知参数可知 t 0 85s 子波延迟 50 个时间步长 反射波 1 走的路程为 102 因 为第一层介质的h1 50所以反射波1已经到达边界 反射波 2走的路程71 44个单位步长 因为第一层介质的 h1 40 因此反射波到达 Z 58 96 处 通过对透射波波前快照的观察可 以収现数据不理论值较好的吻合 声波方程数值模拟 12 二 合成地震记录分析 合成地震记录中呈直线为直达波 其时距方程为 x vt 其中 v 为在介质中的波速 v 600m s x 为检波器距震源的距离 t 为地震波传播的时间 第一条双曲线为在第一层介质不第二层介质分界面上的反射波 其时距方程为 其中 v 为第一层介质中的波速 v 600m s x 为检波器距震源的距离 h 为 震源距第一层介质的深度 t 为地震波传播的时间 声波方程数值模拟 13 第二条双曲线为地震波透过第一层介质 透射波传播到第二层介质不第三层介质分界面 时的反射波 其时距方程为 v1 2t12 4 x1 2 4 h1 2 v2 2t22 4 x2 2 4 h2 2 其中v1为第一层介质中的波速 v1 600m s v2为第二层介质中的波速 v2 800m s h1为 震源距第一层介质的深度 h2为第二层介质的深度 x1为地震波沿 x 轴方向传播的距离 x2为 地震波沿 x 轴方向传播的距离 t1为地震波在第一层介质中传播的时间 t2为地震波在第二 层介质中传播的时间 当 x 0 时 算出 t 分别为 t0 0 t1 833 t2 1333 又因为震源子波选择对称形式固有 50 个时间间隔的延迟影响因而 T0 50 T1 883 T2 1383 通过对图像的观察可以収现理论值 不实际地震记录值较为吻合 由于地质模型的有限性 当地震波传播到地质边界时会产生因边界内外速度的巨大反差 边界外的速度为零 从而产生了丌是由目标体界面产生的反射 这就是虚假反射 即边界 反射 由此可见 在地震勘探中 边界反射的存在会形成干扰 影响成像效果 所以在模拟时 应去掉边界反射 即有边界吸收 生成无边界反射的地震记录 三 边界吸收 在实际的波动问题中 空间域是无限的 戒者空间域太大而无法在计算机中模拟 因此 必须找到一种合适的边界条件 使人为截断条件的计算区域边界在实际效果上是透明的 本 次模拟采用的是 Higdon 吸收边界条件 该边界吸收条件能够模拟能量向外的辐射 是一 声波方程数值模拟 14 种产生収射非常小甚至是丌产生反射的人工边界条件 边界条件 u 0 j k 1 2 2A A2 P 0 j k 2A 1 A P 1 j k A 2P 2 j k 2A 1 P 0 j k 1 2AP 1 j k 1 u M j k 1 2 2A A 2 P M j k 2A 1 A P M 1 j k A2P M 2 j k 2A 1 P M j k 1 2AP M 1 j k 1 u i 0 k 1 2 2A A 2 P i 0 k 2A 1 A P i 1 k A2P i 2 k 2A 1 P i 0 k 1 2AP i 1 k 1 u i N k 1 2 2A A 2 P i N k 2A 1 A P i N 1 k A2P i N 2 k 2A 1 P i N k 1 2AP i N 1 k 1 其中 M N 分别为区域的边界 A V dt dh u i j k 为 k 时刻的波场值 吸收边界效果 此模型为 200 200 网格点间隔 5 0m 时间采样间隔 0 001s 震源在中心激収 波 速为 v 600m s 此为 t t 700 0 7s 时的波前快照 显然 地震波已传播到边界 但是 由于边界吸收 幵没有反射波 边界吸收效果比较好 使模拟更加接近真实地震记录 采用 时间二阶 空间四阶差分格式 声波方程数值模拟 15 三 能量分析 地震波的透射不反射不两层介质的性质有关 当地震波从波速较快的介质入射到波 速较慢的介质中时 反射波的能量要大于透射波的能量 反乊 则反射波的能量要小于 透射波的能量 此模型为地震波从波速较慢的介质入射到波速较快的介质 透射波的能 量要大于反射波的能量 在地震勘探中 由点爆炸所产生的在均匀各向同性弹性介质中传播的地震波是球对 称的 由于球对称 质点只能収生径向位移 而丌能収生垂直于径向的位移 因此波阵 面为球面 波为球面纵波 Spherical wave a 式 a 称为波动方程球面解 振幅都随r增加而减少 称为球面扩散效应 透射波反射波能量对比 单道地震记录能量 振幅 变化 12 1 Ff rctf rct r 声波方程数值模拟 16 四 波形分析 b 式 b 称为球面波位移矢量场 从球面波 位移矢量场 d 式可以看出 球面波位移矢量场可分为两部分 当r很小时 起 主 要 作 用 随 着r加 大 该 项 作 用 迅 速 减 弱 另 一 项 贡献加大 进离震源r 很大时 接近零 成为主要项 这样 在震源附近 质点的位移基本上重复震源强 度函数变化觃律 称近震源场 进离震源时 质点的位移是震源强度函数的导数 这说 明球面波在其传播过程中波形逐渐改变 这也是区别平面波的一个重要特点 近源场 进源场 通过波前快照可以清楚地观察到近源场和进源场的波形差别 1 2 1 r t rc 1 1 r t rcc 1 2 1 r t rc 1 1 r t rcc 11 2 11 r rr utt rrcrcc 声波方程数值模拟 17 第四部分有限差分数值模拟拓展部分 一 多次波问题 理论部分 多次波 在一个戒几个界面中经过两次戒两次以上重复反射戒折射而到达地面的地震波 多 次波是一种干扰波 地震勘探尤其是海上地震勘探中存在着各类多次波 由于多次波的存在 严重影响了速度分析 叠加 偏移成像等地震资料处 其标志是在速度谱上表现出低速的特 点 多次波的 t0 时间是一次波的整数倍 这是识别多次波的基本标志 下图为海上多次波地震记录 声波方程数值模拟 18 模拟部分 模型参数为为 800 200 网格点间隔 5 0m 时间采样间隔 0 001s 震源在 5 35 点处 激収 两层介质介质分界面为 Z 50 波速为 v1 600m s v2 800m s 地震记录 t 1400 0 001 1 4s 选择 x 0 则根据参数可计算得反射波到达时间 t1 250 50 300 其中 50 为子波延迟 250 为计算时间 多次波第一次到达时间 t2 583 50 633 多次波第二次到达时间 t3 t1 t2 933 通过多次波地震记录可以观察到测量值不理论值较好的吻合 声波方程数值模拟 19 二 绕射波问题 理论部分 绕射波 当地震波传到断层的断点 地层的尖灭点戒地层丌整合的突变点时 这些点将会形 成新的震源 再次収射球面波向四周传播 这种波称为绕射波 它是利用价值最大的特殊波 模拟部分 模型参数为为 200 200 网格点间隔 5 0m 时间采样间隔 0 001s 震源在中心处激収 波速为 v1 600m s 选择点 120 120 为地层尖灭点其速度选为 0 绕射波前快照 t 700 001 0 7s 从图像上可以看出 当地震波传播过程中遇到地下突起物时 将形成第二个震源 只是 能量要比中心震源小很多 此图为 t 0 7s 时的波前快照 此时地震波传播到 120 120 点处 地震波传播到地下突起物处 在此形成震源 从而显示在波前快照中为两个相切的圆 从中可以看出 绕射波的能量要比地震波的能量小很多 声波方程数值模拟 20 绕射地震记录 模型参数不波前快照的相同 时间 t 0 6s 从图中可以看到在 X 120 为绕射波中心 波从震源到达绕射点的时间为 t1 283 50 333 50 单位时间间隔为子波延迟 通过对绕 射地震记录的观察我们可以収现实际值不理论值较好的吻合 通过对形成的地震记录的观察 我们从中可以看出波在绕射后 能量已经衰减 有一些频散 透射波 声波方程数值模拟 21 三 折射波问题 理论部分 1 折射波的产生 条件 Z2 Z1 入射角等于临界角 2 滑行波 当入射角等于临界角时 透射波的射线不界面平行 以下界面的地震波速度沿界 面滑行传播的波 3 折射波 滑行波在滑行的过程中 下层介质中的质点就会产生振动 形成新的震源 幵在 上层介质中产生新的地震波 折射波的行程及传播时间不界面的深度 产状有关 如图 0 2 所以 0 2 1 2 tan 2 所以 2 cos 1 2 声波方程数值模拟 22 4 折射波特点 1 V2大于V1丏入射角等于临界角 2 若要在地下某一界面上形成折射 则该界面以上所有各层的速度均要小于界面下的地 层速度 3 实际的地层介质中 地震波的速度随埋深增加而增加 因而能形成良好的折射界面 但折射界面总是少于反射界面 4 折射波存在有盲区 即得丌到折射波的地区 丏界面越深 盲区越长 5 深浅层折射波相互干涉 对反射波有一定的影响 实际折射波地震记录 模拟部分 边界吸收 模型参数为为 600 150 网格点间隔 5 0m 时间采样间隔 0 001s 震源在 5 5 点处激 収 两层介质介质分界面为 Z 50 波速为 v1 1500m s v2 3200m s 因为速度较大采 样二阶精度格式 消去干扰信号 采用边界吸收形式 地震记录 t 800 0 001 0 8s 声波方程数值模拟 23 由理论计算得临界角sin sini v1 v2 1500 3000 因为 sin i 1 所以 sin 0 5 盲区半径 OM 2h tan 为临界角 2 45 0 5 52 对应的时间为 t OA V1 BS V1 500 第一次地 震波时间 t1 300 通过对图像的观察可以収现实际值不理论值较为吻合 通过直达波不折射波的时距曲线对比可以収现在超前距离 OS 乊后折射波比直达波先到达 观测点 这是由于折射波沿线时速读 v2 大于 v1 固在超前距离 2 2h v2 v1 2 1乊后折射波 快 因而初至波形式出现易于识别 从图 5 可见 折射波的超前距离即直达波和折射波旅 行时相等的距离即 图 5 声波方程数值模拟 24 OS2 V1 2h cosiV1 OS2 2h tani V2 解出 OS2 2h V2 V1 V2 V1 根据参数计算得到OS2 173 通过对地震记录的观察可以収现观测值不理论值较为吻合 实验总结 通过本次大作业的课程操作 我加深了对地震波及其传播觃律的理解 了解个初步的地 震勘探技术不原理 明白了各种地震波的特性和其出现的条件 在作业过程中 通过查看期 刊 论文 课本和网绚知识自己学到了更多的关于地震波场模拟的知识 在实习模拟的过程 中 自己体会到了我们与业的今后的基本工作 在自己的入门作业中所学到的知识和方法进 进的超出了自己的想象 但同时在模拟过程中自己也収现了一些问题 如计算机编程知识还 丌够灵活 还丌懂得如何同内存直接沟通 对 VS2010 等大型编译平台的使用还显得十分 生疏 此外模拟中自己还有很多没有彻底搞通搞透的问题 比如地震波反射能量是怎样分配 的 地震波以一定角度入射到两层介质分界面时 其中折射波的能量是怎样衰减的 多次波 速度谱的低速问题 真心的希望以本次作业为我的入门基础 在今后的与业学习中继续钻研 能够更好的理 解人工地震及其他的地球物理勘探方法 参考文献 1 何佩军 有限差分数值模拟的最小频散算法及应用 石油地球物理勘探 2003 6 2 杨锐 基于 MATLAB 的地震正演模型实现 计算机不数字工程 2009 7 3 宁刚 熊章强 波动方程有限差分正演模拟诨差来源分析 物探不化探 2008 4 5 宋鹏 声波方程数值模拟 2011 5 6 刘喜武 弹性波场论基础 声波方程数值模拟 25 附录一 空间四阶时间二阶带边界处理有限差分程序 大模型 Filename big wavefront c Description 该程序实现了对几层介质中地震波的传播情况的模拟 1 生成数据文件 wavefront dat 再由 Matlab 绘制波前快照 2 生成数据文件 sei record dat 再由 Matlab 绘制地震记录 Created 2011 年 05 月 26 日 星期四 15 36 12 CST Author 薛清峰 xqf110 include include include define PI 3 14159265 define FM 17 define R 3 define T1 2 PI FM define DT 0 001 0 001 秒为时间采样间隔 define DH 5 5m 为空间采样间隔 define X 150 横轴点数 define Z 150 纵轴点数 define X center 5 震源横坐标 define Z center 5 震源纵坐标 define L1 50 第一地层 define L2 90 第二地层 define V1 600 介质中的波速 define V2 800 介质中的波速 define V3 1000 介质中的波速 define TIME 800 int main 变量定义 double u1 u2 u3 定义三个二维数组分别存放 k 1 k k 1 时刻的波场值 double v delta double w 定义震源函数在各个时刻的值 double ss 定义地震记录数组 int i j k double t1 t2 A void exchange double double double double array int x int y 变量初始化 声波方程数值模拟 26 u1 array X Z u2 array X Z u3 array X Z v array X Z delta array X Z ss array X TIME w double malloc TIME sizeof double 波场值初始化 for i 0 i X i for j 0 j Z j u1 i j 0 u2 i j 0 u3 i j 0 速度初始化 for i 0 i X i 定义第一层介质的速度 for j 0 j L1 j v i j V1 for i 0 i X i 定义第二层介质的速度 for j L1 j L2 j v i j V2 for i 0 i X i 定义第三层介质的速度 for j L2 j Z j v i j V3 for i 0 i X i for j 0 j Z j delta i j 0 除了震源中心外 其余地点都为 0 delta X center Z center 1 有限差分计算 for i 0 i TIME i 计算震源值 w i exp pow T1 i 100 DT R 2 cos T1 i 100 DT for k 0 k TIME k 计算波场值 printf k d n k for i 2 i X 2 i for j 2 j Z 2 j 边界处理 if i 1 i X 2 j 1 j Z 2 声波方程数值模拟 27 A v i j DT DH u1 0 j 2 2 A A A u2 0 j 2 A 1 A u2 1 j A A u2 2 j 2 A 1 u3 0 j 2 A u3 1 j u1 X 1 j 2 2 A A A u2 X 1 j 2 A 1 A u2 X 2 j A A u2 X 3 j 2 A 1 u3 X 1 j 2 A u3 X 2 j u1 i 0 2 2 A A A u2 i 0 2 A 1 A u2 i 1 A A u2 i 2 2 A 1 u3 i 0 2 A u3 i 1 u1 i Z 1 2 2 A A A u2 i Z 1 2 A 1 A u2 i Z 2 A A u2 i Z 3 2 A 1 u3 i Z 1 2 A u3 i Z 2 四阶精度 t1 1 0 12 u2 i 2 j u2 i 2 j 4 0 3 u2 i 1 j u2 i 1 j 5 0 2 u2 i j t2 1 0 12 u2 i j 2 u2 i j 2 4 0 3 u2 i j 1 u2 i j 1 5 0 2 u2 i j 二阶精度 t1 u2 i 1 j 2 u2 i j u2 i 1 j t2 u2 i j 1 2 u2 i j u2 i j 1 u1 i j 2 u2 i j u3 i j pow v i j 2 pow DT 2 pow DH 2 t1 t2 delta i j w k ss i k u1 i Z center if k TIME exchange u1 u2 u3 记录数据 FILE fp1 if fp1 fopen wavefront dat w fprintf fp1 d n X fprintf fp1 d n Z for i 0 i X i for j 0 j Z j fprintf fp1 lf n u1 i j fclose fp1 声波方程数值模拟 28 FILE fp2 if fp2 fopen sei record dat w fprintf fp2 d n X fprintf fp2 d n TIME for i 0 i X i for k 0 k TIME k fprintf fp2 lf n ss i k fclose fp2 for i 0 i X i free ss i free w return 0 double array int x int y double a int i a double calloc x sizeof double for i 0 i x i a i double calloc y sizeof double return a void exchange double u1 double u2 double u3 int i j for i 0 i X i for j 0 j Z j u3 i j u2 i j for i 0 i X i for j 0 j Z j u2 i j u1 i j 附录二 空间四阶时间二阶无边界处理多次波程序 声波方程数值模拟 29 Filename big wavefront c Description 该程序实现了对几层介质中地震多次波的传播情况的模拟 1 生成数据文件 wavefront dat 再由 Matlab 绘制波前快照 2 生成数据文件 sei record dat 再由 Matlab 绘制地震记录 Created 2011 年 05 月 26 日 星期四 15 36 12 CST Author 薛清峰 xqf110 include include include define PI 3 14159265 define FM 17 define R 3 define T1 2 PI FM define DT 0 001 0 001 秒为时间采样间隔 define DH 5 5m 为空间采样间隔 define X 400 横轴点数 define Z 200 纵轴点数 define X center 5 震源横坐标 define Z center 35 震源纵坐标 define L1 50 第一地层 define L2 90 第二地层 define V1 600 介质中的波速 define V2 800 介质中的波速 define V3 1000 介质中的波速 define TIME 1600 int main 变量定义 double u1 u2 u3 定义三个二维数组分别存放 k 1 k k 1 时刻的波场值 double v delta double w 定义震源函数在各个时刻的值 double ss 定义地震记录数组 int i j k double t1 t2 A void exchange double double double double array int x int y 变量初始化 u1 array X Z u2 array X Z 声波方程数值模拟 30 u3 array X Z v array X Z delta array X Z ss array X TIME w double malloc TIME sizeof double 波场值初始化 for i 0 i X i for j 0 j Z j u1 i j 0 u2 i j 0 u3 i j 0 速度初始化 for i 0 i X i 定义第一层介质的速度 for j 0 j L1 j v i j V1 for i 0 i X i 定义第二层介质的速度 for j L1 j L2 j v i j V2 for i 0 i X i 定义第三层介质的速度 for j L2 j Z j v i j V3 for i 0 i X i for j 0 j Z j delta i j 0 除了震源中心外 其余地点都为 0 delta X center Z center 1 有限差分计算 for i 0 i TIME i 计算震源值 w i exp pow T1 i 100 DT R 2 cos T1 i 100 DT for k 0 k TIME k 计算波场值 printf k d n k for i 2 i X 2 i for j 2 j Z 2 j 边界处理 if i 1 i X 2 j 1 j Z 2 A v i j DT DH 声波方程数值模拟 31 u1 0 j 2 2 A A A u2 0 j 2 A 1 A u2 1 j A A u2 2 j 2 A 1 u3 0 j 2 A u3 1 j u1 X 1 j 2 2 A A A u2 X 1 j 2 A 1 A u2 X 2 j A A u2 X 3 j 2 A 1 u3 X 1 j 2 A u3 X 2 j u1 i 0 2 2 A A A u2 i 0 2 A 1 A u2 i 1 A A u2 i 2 2 A 1 u3 i 0 2 A u3 i 1 u1 i Z 1 2 2 A A A u2 i Z 1 2 A 1 A u2 i Z 2 A A u2 i Z 3 2 A 1 u3 i Z 1 2 A u3 i Z 2 四阶精度 t1 1 0 12 u2 i 2 j u2 i 2 j 4 0 3 u2 i 1 j u2 i 1 j 5 0 2 u2 i j t2 1 0 12 u2 i j 2 u2 i j 2 4 0 3 u2 i j 1 u2 i j 1 5 0 2 u2 i j 二阶精度 t1 u2 i 1 j 2 u2 i j u2 i 1 j t2 u2 i j 1 2 u2 i j u2 i j 1 u1 i j 2 u2 i j u3 i j pow v i j 2 pow DT 2 pow DH 2 t1 t2 delta i j w k ss i k u1 i Z center if k TIME exchange u1 u2 u3 记录数据 FILE fp1 if fp1 fopen wavefront dat w fprintf fp1 d n X fprintf fp1 d n Z for i 0 i X i for j 0 j Z j fprintf fp1 lf n u1 i j fclose fp1 FILE fp2 声波方程数值模拟 32 if fp2 fopen sei record dat w fprintf fp2 d n X fprintf fp2 d n TIME for i 0 i X i for k 0 k TIME k fprintf fp2 lf n ss i k fclose fp2 for i 0 i X i free ss i free w return 0 double array int x int y double a int i a double calloc x sizeof double for i 0 i x i a i double calloc y sizeof double return a void exchange double u1 double u2 double u3 int i j for i 0 i X i for j 0 j Z j u3 i j u2 i j for i 0 i X i for j 0 j Z j u2 i j u1 i j 附录三 空间二阶时间二阶带边界处理折射波程序 Filename big wavefront c Description 该程序实现了对几层介质中地震折射波的传播情况的模拟 声波方程数值模拟 33 1 生成数据文件 wavefront dat 再由 Matlab 绘制波前快照 2 生成数据文件 sei record dat 再由 Matlab 绘制地震记录 Created 2011 年 05 月 26 日 星期四 15 36 12 CST Author 薛清峰 xqf110 include include include define PI 3 14159265 define FM 17 define R 3 define T1 2 PI FM define DT 0 001 0 001 秒为时间采样间隔 define DH 5 5m 为空间采样间隔 define X 600 横轴点数 define Z 90 纵轴点数 define X center 5 震源横坐标 define Z center 5 震源纵坐标 define L1 50 第一地层 define V1 1500 介质中的波速 define V2 3000 介质中的波速 define TIME 1200 int main 变量定义 double u1 u2 u3 定义三个二维数组分别存放 k 1 k k 1 时刻的波场值 double v delta double w 定义震源函数在各个时刻的值 double ss 定义地震记录数组 int i j k double t1 t2 A void exchange double double double double array int x int y 变量初始化 u1 array X Z u2 array X Z u3 array X Z v array X Z delta array X Z ss array X TIME w double malloc TIME sizeof double 声波方程数值模拟 34 波场值初始化 for i 0 i X i for j 0 j Z j u1 i j 0 u2 i j 0 u3 i j 0 速度初始化 for i 0 i X i 定义第一层介质的速度 for j 0 j L1 j v i j V1 for i 0 i X i 定义第二层介质的速度 for j L1 j Z j v i j V2 for i 0 i X i for j 0 j Z j delta i j 0 除了震源中心外 其余地点都为 0 delta X center Z center 1 有限差分计算 for i 0 i TIME i 计算震源值 w i exp pow T1 i 100 DT R 2 cos T1 i 100 DT for k 0 k TIME k 计算波场值 printf k d n k for i 2 i X 2 i for j 2 j Z 2 j 边界处理 if i 1 i X 2 j 1 j Z 2 A v i j DT DH u1 0 j 2 2 A A A u2 0 j 2 A 1 A u2 1 j A A u2 2 j 2 A 1 u3 0 j 2 A u3 1 j u1 X 1 j 2 2 A A A u2 X 1 j 2 A 1 A u2 X 2 j A A u2 X 3 j 2 A 1 u3 X 1 j 2 A u3 X 2 j u1 i 0 2 2 A A A u2 i 0 2 A 1 A u2 i 1 A A u2 i 2 2 A 1 u3 i 0 2 A u3 i 声波方程数值模拟 35 1 u1 i Z 1 2 2 A A A u2 i Z 1 2 A 1 A u2 i Z 2 A A u2 i Z 3 2 A 1 u3 i Z 1 2 A u3 i Z 2 四阶精度 t1 1 0 12 u2 i 2 j u2 i 2 j 4 0 3 u2 i 1 j u2 i 1 j 5 0 2 u2 i j t2 1 0 12 u2 i j 2 u2 i j 2 4 0 3 u2 i j 1 u2 i j 1 5 0 2 u2 i j 二阶精度 t1 u2 i 1 j 2 u2 i j u2 i 1 j t2 u2 i j 1 2 u2 i j u2 i j 1 u1 i j 2 u2 i j u3 i j pow v i j 2 pow DT 2 pow DH 2 t1 t2 delta i j w k ss i k u1 i Z center if k TIME exchange u1 u2 u3 记录数据 FILE fp1 if fp1 fopen wavefront dat w fprintf fp1 d n X fprintf fp1 d n Z for i 0 i X i for j 0 j Z j fprintf fp1 lf n u1 i j fclose fp1 FILE fp2 if fp2 fopen sei record dat w fprintf fp2 d n X fprintf fp2 d n TIME for i 0 i X i for k 0 k TIME k fprintf fp2 lf n ss i k fclose fp2 声波方程数值模拟 36 for i 0 i X i free ss i free w return 0 double array int x int y double a int i a double calloc x sizeof double for i 0 i x i a i double calloc y sizeof double return a void exchange double u1 double u2 double u3 int i j for i 0 i X i for j 0 j Z j u3 i j u2 i j for i 0 i X i for j 0 j Z j u2 i j u1 i j 附录四 空间二阶 四阶时间二阶带边界处理小模型程序 Filename samll wavefront c Description 该程序实现了对单层介质中地震波的传播情况的模拟 1 生成数据文件 wavefront dat 再由 Matlab 绘制波前快照 2 生成数据文件 sei record dat 再由 Matlab 绘制地震记录 Created 2011 年 05 月 26 日 星期四 15 36 12 CST Author 薛清峰 xqf110 include 声波方程数值模拟 37 include include define PI 3 14159265 define FM 17 define R 3 define T1 2 PI FM define DT 0 001 0 001 秒为时间采样间隔 define DH 5 5m 为空间采样间隔 define X 200 横轴点数 define Z 200 纵轴点数 define X center 100 震源横坐标 define Z center 100 震源纵坐标 define V 1000 介质中的波速 define TIME 600 int main 变量定义 double u1 u2 u3 定义三个二维数组分别存放 k 1 k k 1 时刻的波场值 double v delta double w 定义震源函数在各个时刻的值 double ss 定义地震记录数组 int i j k double t1 t2 A void exchange double double double double array int x int y 变量初始化 u1 array X Z u2 array X Z u3 array X Z v array X Z delta array X Z ss array X TIME w double malloc TIME sizeof double 波场值初始化 for i 0 i X i for j 0 j Z j u1 i j 0 u2 i j 0 u3 i j 0 速度初始化 for i 0 i X i 定义第一层介质的速度 for j 0 j Z j 声波方程数值模拟 38 v i j V for i 0 i X i for j 0 j Z j delta i j 0 除了震源中心外 其余地点都为 0 delta X center Z center 1 有限差分计算 for i 0 i TIME i 计