地球物理计算方法课件：绪论2
合集下载
计算地球物理学
计算地球物理学计算地球物理学是地球科学中的一个重要分支，它研究地球内部的物质组成、结构及其运动规律。
通过计算地球物理学的方法，我们可以揭示地球内部的奥秘，为地质灾害预测、资源勘探和环境保护等领域提供科学依据。
地球物理学的计算方法主要包括地震学、重力学和磁力学等。
地震学是通过地震波的传播和接收来研究地球内部结构的一门学科。
通过分析地震波的传播速度和路径，地震学家可以推断出地球内部的物质性质和分布情况。
重力学则是通过测量地球重力场的变化来研究地球内部的密度分布。
磁力学则是通过测量地球磁场的变化来研究地球内部的磁性物质分布。
地震学中的计算方法主要包括地震波传播速度的计算和地震震源的定位。
地震波传播速度的计算是地震学研究的重要内容之一，它可以揭示地球内部的物质性质和结构。
地震波传播速度的计算方法主要有声波和弹性波传播理论，通过求解波动方程可以得到地震波传播速度的分布。
地震震源的定位是通过测量地震波到达不同测点的时间差来确定地震震源的位置，而计算地震震源的定位则需要运用到三角测量和时间差测量等方法。
重力学中的计算方法主要包括重力场的测量和重力异常的计算。
重力场的测量可以通过使用重力仪器进行，通过测量重力加速度的变化可以得到重力场的分布。
而重力异常的计算则是通过将实测重力场与正常重力场之间的差异进行计算，从而得到地下物质的密度分布情况。
磁力学中的计算方法主要包括磁场的测量和磁异常的计算。
磁场的测量可以通过使用磁力仪器进行，通过测量磁场的变化可以得到磁场的分布情况。
而磁异常的计算则是通过将实测磁场与正常磁场之间的差异进行计算，从而得到地下磁性物质的分布情况。
计算地球物理学在地球科学研究和实践中发挥着重要的作用。
通过计算地球物理学的方法，我们可以了解地球内部的物质组成和结构，揭示地球演化的历史和过程。
同时，计算地球物理学还可以应用于地质灾害的预测和监测，通过分析地震活动和重力、磁场的变化，可以提前预警地震、火山喷发等灾害事件。
地球物理勘探概论普通物探方法PPT课件
第7页/共109页
重力勘探是在地表观测由于地下岩矿石存在密度差异 而引起的重力场的变化，通过对观测资料的处理和解 释，达到反求地下构造的一种勘探方法。
该勘探方法勘探的地球物理前提条件是：地下岩矿石 存在密度差异。
所依据讨论的物理场是重力场。重力场是一种位场。 即场的变化只与空间变化有关，与时间无关。
第2页/共109页
第2章 普通物探方法
大地电磁勘探技术是比较理想的方法之一，并且随着大地电磁仪器的发展，基础理论的 进步以及处理方法的改进，大地电磁方法日臻完善。因此，20世纪90年代以来，大地 电磁方法在国内外都得到了广泛的应用。
➢ 重、磁、电、地震联合反演 重磁约束反演是近年来重磁应用的新成果。其基本原理是利用几条地震及电法剖面做约
国内外技术现状：目前国内重力勘探队伍主要集中 于地矿部门，石油部门的东方地球物理勘探公司有完 备的集野外采集、处理及解释为一体的专业队伍。从 在胜利油田施工的队伍－北京勘察技术工程公司来看， 其测量精度达到微伽级，代表了国内领先水平。从国 外研究看，重力勘探除应用于盆地及盆地深层的勘探 外，井中重力测井也取得了一定进展。
第23页/共109页
第2章 普通物探方法-§2-1 重力勘探
陕甘宁盆地桌子山段布格重力异常与不同高度延拓的 重力异常图，反映了本段深浅层构造方向不相同，不 是同一期构造运动的产物。
第20页/共109页
第2章 普通物探方法-§2-1 重力勘探
（2） 研究沉积盖层内部构造 绝大多数的油气藏都在沉积盖层中，多数沉积盖层构
造与基底构造、断裂活动、乃至火成岩活动有关，因此 通过研究结晶基岩可以间接寻找沉积岩构造。
主要的情况有：①盖层构造与基底高点具有继承关系； ②根据断层的存在及其性质可以推断与其相关的构造， 如滚动背斜；③沉积盆地后期的岩浆侵入可以在盖层中 造成正构造形态，此时，沉积盖层构造的核部是火成岩 体。
重力勘探是在地表观测由于地下岩矿石存在密度差异 而引起的重力场的变化，通过对观测资料的处理和解 释，达到反求地下构造的一种勘探方法。
该勘探方法勘探的地球物理前提条件是：地下岩矿石 存在密度差异。
所依据讨论的物理场是重力场。重力场是一种位场。 即场的变化只与空间变化有关，与时间无关。
第2页/共109页
第2章 普通物探方法
大地电磁勘探技术是比较理想的方法之一，并且随着大地电磁仪器的发展，基础理论的 进步以及处理方法的改进，大地电磁方法日臻完善。因此，20世纪90年代以来，大地 电磁方法在国内外都得到了广泛的应用。
➢ 重、磁、电、地震联合反演 重磁约束反演是近年来重磁应用的新成果。其基本原理是利用几条地震及电法剖面做约
国内外技术现状：目前国内重力勘探队伍主要集中 于地矿部门，石油部门的东方地球物理勘探公司有完 备的集野外采集、处理及解释为一体的专业队伍。从 在胜利油田施工的队伍－北京勘察技术工程公司来看， 其测量精度达到微伽级，代表了国内领先水平。从国 外研究看，重力勘探除应用于盆地及盆地深层的勘探 外，井中重力测井也取得了一定进展。
第23页/共109页
第2章 普通物探方法-§2-1 重力勘探
陕甘宁盆地桌子山段布格重力异常与不同高度延拓的 重力异常图，反映了本段深浅层构造方向不相同，不 是同一期构造运动的产物。
第20页/共109页
第2章 普通物探方法-§2-1 重力勘探
（2） 研究沉积盖层内部构造 绝大多数的油气藏都在沉积盖层中，多数沉积盖层构
造与基底构造、断裂活动、乃至火成岩活动有关，因此 通过研究结晶基岩可以间接寻找沉积岩构造。
主要的情况有：①盖层构造与基底高点具有继承关系； ②根据断层的存在及其性质可以推断与其相关的构造， 如滚动背斜；③沉积盆地后期的岩浆侵入可以在盖层中 造成正构造形态，此时，沉积盖层构造的核部是火成岩 体。
地球物理解释基础ppt课件
第19章介绍AVO方法
• 什么是AVO?——研究CMP道集内相对振幅，称作振幅随炮检
距变化的分析（AVO）。研究相对振幅随反射角的变化，这种方法称 作振幅随入射角变化的分析（AVA）
• AVO分析能解决什么地质问题 —碎屑岩气藏直接烃类指示 —在碳酸盐油藏中可能识别孔隙发育带
• 随炮检距变化的反射系数 ——AVO的计算
地球物理解释基础 (2)
第14章解释与盐构造有关的圈闭
• 许多重要油田 和盐圈闭联系在一起——著名的墨西哥湾、美国几
个洲、加拿大、北海、北非、德国 、里海地区 都存在盐圈闭。塔里木盆 地的克拉2气田也和盐圈闭有关
• 盐与众不同 ——有较低的密度和较高的地震速度 • “漂浮”状侵入到沉积物之下——侵入体产生各种盐体形状 ，盐
层析成像技术
• 层析的目标是求解每个面元的速度 • 层析成像方法的第一步是从未叠加的地震资料上或直接从野
外观测值拾取旅行时 • 建立初始模型作射线追踪，由初始模型计算的波至时间与观
测值进行比较 • 根据模型值与观测值之差对模型进行修改 • 拾取旅行时-建立模型-模型正演-测量时间差-修改模型——这
正演模拟和反演的关系
正演模拟——用一个数学关系式，对给
出的一组模型参数合成地下响应。
反演或“反演模拟” ——与正演 模拟
“相反”的过程。对一个给出的数据集， 寻
求定义一个与观测数据相符的地质模型
从数学上讲，反问题由于比方程式更多 未知数的存在，能够引起不确定性，产生 多解，所以反演的多解性是固有的
反演的多解性（非唯一性、不确定性）
1D, 2D, 3D需要将物理模 型材料校正为成比例的模 型
当非均质体的尺度与 Fresnel带相比很大时， 一般是可以应用的。通
地球物理正演方法
值越小，计算值与理论值越接近。
矛盾：减小步长h将成倍增加计算节点数目，增加计算机内存需
求和计算时间。降低了效率，增加了费用
解决计算速度与精度矛盾的较好方法：采用变步长，即在
近区将网格分得密些，远区影响较小，可分得稀些。
弹性波场计算举例 1、反射地震中波传播方程
密度不均匀介质弹性波标量波动方程
激发问题
模拟一个或少数几个小异常体的响应时，该方法比较经济 多用于3D数值模拟
假设大地电磁场的源是来自高空的 垂直入射到地面的平面电磁波，则 频域中无源麦克斯韦方程组
总 场
表示地下任一点处的实际电导 率值，且有
1，异常体外 = 2，异常体内
2(r)
1
V
r表示矢径
定义：一次场为均匀地球场，并以上角标p表示，则一次场也满足无源 麦克斯韦方程组，有
在各向同性均匀介质、平面波入射假设条件下，标量波动方程
传播问题
在二维情况下，
zz|z=0=ux=0，zx|z=0=uz=0
初始条件 （自由表面）边界条件
2、区域离散化
采用正方形网格元进行网格划分， 步长h；m，n为当前网格节点的横 向及垂向编号;l时间取样号
差分方程式
利用差分方程式，由上至下，由左至右并随时标l增加计算 空间任一点(m，n)的波场um，n，l+1便得到波传播图像， um，0，l是地面直达波和反射波场的合成记录。
2、将势场u展成某种简单函数和系数的线性组合 假定，单元内势可用线性(一阶)方程表示，有 V=a+bx+cy 沿三角元边缘势V可以由相应两角点势值线性内插而 来，如果两个三角元共用一条边，则位势在跨单元时 保持连续。 为求各系数，设三个顶点上势为V1，V2，V3
矛盾：减小步长h将成倍增加计算节点数目，增加计算机内存需
求和计算时间。降低了效率，增加了费用
解决计算速度与精度矛盾的较好方法：采用变步长，即在
近区将网格分得密些，远区影响较小，可分得稀些。
弹性波场计算举例 1、反射地震中波传播方程
密度不均匀介质弹性波标量波动方程
激发问题
模拟一个或少数几个小异常体的响应时，该方法比较经济 多用于3D数值模拟
假设大地电磁场的源是来自高空的 垂直入射到地面的平面电磁波，则 频域中无源麦克斯韦方程组
总 场
表示地下任一点处的实际电导 率值，且有
1，异常体外 = 2，异常体内
2(r)
1
V
r表示矢径
定义：一次场为均匀地球场，并以上角标p表示，则一次场也满足无源 麦克斯韦方程组，有
在各向同性均匀介质、平面波入射假设条件下，标量波动方程
传播问题
在二维情况下，
zz|z=0=ux=0，zx|z=0=uz=0
初始条件 （自由表面）边界条件
2、区域离散化
采用正方形网格元进行网格划分， 步长h；m，n为当前网格节点的横 向及垂向编号;l时间取样号
差分方程式
利用差分方程式，由上至下，由左至右并随时标l增加计算 空间任一点(m，n)的波场um，n，l+1便得到波传播图像， um，0，l是地面直达波和反射波场的合成记录。
2、将势场u展成某种简单函数和系数的线性组合 假定，单元内势可用线性(一阶)方程表示，有 V=a+bx+cy 沿三角元边缘势V可以由相应两角点势值线性内插而 来，如果两个三角元共用一条边，则位势在跨单元时 保持连续。 为求各系数，设三个顶点上势为V1，V2，V3
计算物理基础课件
第二十页，共609页。
1.2 计算物理的起源、形成与发展
1963年，美国的Beini，Alder等人开始编辑出版《计 算物理方法》丛书，内容涉及统计物理、量子力学、流 体力学、核物理、天体物理、固体物理、等离子体物理、 地球物理和大气环流等。
1966年，Journal of Computational Physics在美国创刊；1969年，Computer Physics Communication在西欧创刊。
多人认为是计算物理的正式起点，因为它提出了 许多问题，带来了当时谁也未曾想到的重大发 展。
第十八页，共609页。
1.2 计算物理的起源、形成与发展
从此，物理问题的计算与计算机相互促进，开 始蓬勃发展。
1950年，全世界还只有15台计算机，
到1962年9月，仅美国就有了16817台。
现在的计算机不计其数！
理的问题），并于1952年夏天设计了一个计算机实验， 一年后，在当时用来进行氢弹设计的MANIAC计算机上 实现。
第十七页，共609页。
1.2 计算物理的起源、形成与发展
1954年11月，费米逝世，他的合作者继续工作， 于 1955 年 5 月 写 出 Los Alamos 研 究 报 告 LA1940。这篇秘密报告历经多年、解密后被正式收 入《费米全集》。这篇具有重大意义的报告，被许
第八页，共609页。
1.1 什么是计算物理？
➢ 计算物理是伴随着电子计算机的出现和发 展而逐步形成的一门新兴的边缘学科。
➢ 是以电子计算机为工具、采用数学方法解决 物理问题的应用科学。
➢ 是物理、数学和计算机三者相结合的产物。
第九页，共609页。
1.1 什么是计算物理？
计算物理中的“计算”，不是上物理课做习 题时进行的那种简单计算；不是用古典的数学物理 方法来完成的计算；
1.2 计算物理的起源、形成与发展
1963年，美国的Beini，Alder等人开始编辑出版《计 算物理方法》丛书，内容涉及统计物理、量子力学、流 体力学、核物理、天体物理、固体物理、等离子体物理、 地球物理和大气环流等。
1966年，Journal of Computational Physics在美国创刊；1969年，Computer Physics Communication在西欧创刊。
多人认为是计算物理的正式起点，因为它提出了 许多问题，带来了当时谁也未曾想到的重大发 展。
第十八页，共609页。
1.2 计算物理的起源、形成与发展
从此，物理问题的计算与计算机相互促进，开 始蓬勃发展。
1950年，全世界还只有15台计算机，
到1962年9月，仅美国就有了16817台。
现在的计算机不计其数！
理的问题），并于1952年夏天设计了一个计算机实验， 一年后，在当时用来进行氢弹设计的MANIAC计算机上 实现。
第十七页，共609页。
1.2 计算物理的起源、形成与发展
1954年11月，费米逝世，他的合作者继续工作， 于 1955 年 5 月 写 出 Los Alamos 研 究 报 告 LA1940。这篇秘密报告历经多年、解密后被正式收 入《费米全集》。这篇具有重大意义的报告，被许
第八页，共609页。
1.1 什么是计算物理？
➢ 计算物理是伴随着电子计算机的出现和发 展而逐步形成的一门新兴的边缘学科。
➢ 是以电子计算机为工具、采用数学方法解决 物理问题的应用科学。
➢ 是物理、数学和计算机三者相结合的产物。
第九页，共609页。
1.1 什么是计算物理？
计算物理中的“计算”，不是上物理课做习 题时进行的那种简单计算；不是用古典的数学物理 方法来完成的计算；
地球物理方法
地球物理方法:利用各种仪器在地表观测地壳上的各种物理现象，从而推断、了解地下的地质构造特点，寻找可能的储油构造。
它是一种间接找油的方法。
特点：精度和成本均高于地质法，但低于钻探方法。
地震勘探：就是利用人工方法激发的地震波(弹性波)，研究地震波在地层中传播的规律，以查明地下的地质构造，从而来确定矿藏（包括油气、矿石、水、地热资源等）等的位置，以及获得工程地质信息。
地震波运动学：研究在地震波传播过程中的地震波波前的空间位置与其传播时间的关系，即研究波的传播规律，以及这种时空关系与地下地质构造的关系。
地震波动力学：研究地震波在传播过程中波形、振幅、频率、相位等特征的及其变化规律，以及这些变化规律与地下的地层结构，岩石性质及流体性质之间存在的联系。
地震波：是一种在岩层中传播的，频率较低(与天然地震的频率相近)的波，弹性波在岩层中传播的一种通俗说法。
地震波由一个震源激发。
地震子波：爆炸产生的是一个延续时间很短的尖脉冲，这一尖脉冲造成破坏圈、塑性带，最后使离震源较远的介质产生弹性形变，形成地震波，地震波向外传播一定距离后，波形逐渐稳定，成为一个具有2-3个相位（极值）、延续时间60-100毫秒的地震波，称为地震子波。
地震子波看作组成一道地震记录的基本元素。
波前：振动刚开始与静止时的分界面，即刚要开始振动的那一时刻。
射线：是用来描述波的传播路线的一种表示。
在一定条件下，认为波及其能量是沿着一条“路径”从波源传到所观测的一点P。
这是一条假想的路径，也叫波线。
射线总是与波阵面垂直，波动经过每一点都可以设想有这么一条波线。
振动图和波剖面：某点振动随时间的变化的曲线称为振动曲线，也称振动图。
地震勘探中，沿测线画出的波形曲线，也称波剖面。
折射波：当入射波大于临界角时，出现滑行波和全反射。
在分界面上的滑行波有另一种特性，即会影响第一界面，并激发新的波。
在地震勘探中，由滑行波引起的波叫折射波，也叫做首波。
入射波以临界角或大于临界角入射高速介质所产生的波滑行波：由透射定律可知，如果V2>V1 ，即sinθ2 > sinθ1 ，θ2 > θ1。
地球物理计算方法课件：上机课2
项式，然后将f(x),P5(x),P10(x),P20(x)画在同一张图上进 行对比分析。（P30，例6，龙格现象） 3. 请自己编写一个Newton插值函数，并重复第2题的步骤， 在此基础上对比分析Lagrange与Newton插值的结果。(选作)
第一次上机作业点评：课后作业
1. 调 用 系 统 函 数 对 课 堂 练 习 2 中 的 f(x) 进 行 分 段 线 性 插 值 、 分 段 三 次 Hermite插值及三次样条插值，并与Lagrange/Newton插值的结果进行对 比分析。
Matlab符号积分函数
Matlab数值积分函数
数值积分
梯形求积公式：
Z = trapz(x,y)
x、y为求积函数一系列已知的节点及函数值， Z返回积分的近似值。
Matlab数值积分函数
自适应复化Simpson公式： [I,n]=quad('fname',a,b,tol,trace)
说明： ① 返回参数I即定积分值，n为被积函数的调用次数； ② fname是被积函数名，定义被积函数须用数值运算符； ③ a和b分别是定积分的下限和上限； ④ tol用来控制积分精度，缺省时取tol=10-6； ⑤ trace控制是否展现积分过程，若取非0则展现积分过程，
取0则不展现，缺省时取trace=0。
Matlab数值积分函数
二重积分： I=dblquad(fun,xmin,xmax,ymin,ymax,tol ,trace)
三重积分： I = triplequad(fun,xmin,xmax,ymin,ymax,zmin,zmax,tol)
高斯积分： I = quadgk(fun,a,b)
格式：f=poly2sym(p)
第一次上机作业点评：课后作业
1. 调 用 系 统 函 数 对 课 堂 练 习 2 中 的 f(x) 进 行 分 段 线 性 插 值 、 分 段 三 次 Hermite插值及三次样条插值，并与Lagrange/Newton插值的结果进行对 比分析。
Matlab符号积分函数
Matlab数值积分函数
数值积分
梯形求积公式：
Z = trapz(x,y)
x、y为求积函数一系列已知的节点及函数值， Z返回积分的近似值。
Matlab数值积分函数
自适应复化Simpson公式： [I,n]=quad('fname',a,b,tol,trace)
说明： ① 返回参数I即定积分值，n为被积函数的调用次数； ② fname是被积函数名，定义被积函数须用数值运算符； ③ a和b分别是定积分的下限和上限； ④ tol用来控制积分精度，缺省时取tol=10-6； ⑤ trace控制是否展现积分过程，若取非0则展现积分过程，
取0则不展现，缺省时取trace=0。
Matlab数值积分函数
二重积分： I=dblquad(fun,xmin,xmax,ymin,ymax,tol ,trace)
三重积分： I = triplequad(fun,xmin,xmax,ymin,ymax,zmin,zmax,tol)
高斯积分： I = quadgk(fun,a,b)
格式：f=poly2sym(p)
地球物理计算方法课件：第二章_数值积分 1
平均高度
a
bx
左矩形公式 f ( ) f (a )
b
a f (x)dx (b a) f (a)
y
右矩形公式 f ( ) f (b )
b
a f (x)dx (b a) f (b)
y
y=f(x)
y=f(x)
11
a
bx a
bx
中矩形公式
f ( ) f ( a b)
2
b
ab
a f (x)dx (b a) f ( 2 )
0
3
计算 f (x) 1, x, x 2 , x3 , x 4 , e x 时积分结果并与准确值进行
比较.
解:梯形公式和Simpson公式的计算结果与准确值比较如下表所示
f(x) 准确值 梯形公式 辛卜生公式
1 x x2
x3 x4
ex
2 2 2.67 4 6.40 6.389
22 4
8 16 8.389
b
1dx
a
b
a，右
( A0
A1 )
当f (x) x时，左 b xdx 1 (b2 a2 )， 右 [A0a A1b]；
a
2
得到：
A0
A1
b
2
a；
例 若对于给定的一组求积节点 xk (k 0,1, 2, , n) 相应的求积公式
b
n
f (x)dx
a
Ak f (xk )
k 0
至少具有n次代数精度，试确定其求积系数．
xnn An
bn1 (n
an1 1) !
其系数矩阵
1 1
x0
x1
x02
x12
x0n x1n
第2章 地球物理中常用数值解法的基本原理-1
i 1 n
n
椭圆型—位场方程 抛物型—扩散方程、 热传导方程
u
xi xi
Biuxi Cu F 0
i 1 n n
ux1 x1 uxi xi Biuxi Cu F 0 双曲型—波动方程
i 2 i 1
第一节 偏微分方程的有限差分解法
2.2 椭圆型偏微分方程的有限差分解法 用差分法将连续问题离散化的步骤： （1）首先对求解区域作网格剖分，用有限个网格节点代替连 续区域；
U xi 1, y j U xi , y j U xi 1, y j U xi , y j
U xi 1, y j U xi 1, y j 2U xi , y j
U xi 1 , y j 2U xi , y j U xi 1 , y j h12
U f x, y ,
hUij Ui 1, j 2Uij Ui 1, j h
2 1
x, y
Ui , j 1 2Uij Ui , j 1 h
2 2
fij
2 截断误差为 O h
第一节 偏微分方程的有限差分解法
1.2 椭圆型偏微分方程的有限差分解法 1.2.1 差分格式—— 九点差分格式
2 1 U xi , y j 2 U xi 1 , y j U xi , y j h1 h1 2 x 2! x 3 4 5 U x , y U x , y 1 1 1 U xi , y j 5 i j i j 3 4 6 h h h O h 1 1 1 1 3 4 5 3! x 4! x 5! x 2 1 U xi , y j 2 U xi 1 , y j U xi , y j h1 h1 2 x 2! x 3 4 5 1 U xi , y j 3 1 U xi , y j 4 1 U xi , y j 5 6 h h h O h 1 1 1 1 3! x3 4! x 4 5! x 5
n
椭圆型—位场方程 抛物型—扩散方程、 热传导方程
u
xi xi
Biuxi Cu F 0
i 1 n n
ux1 x1 uxi xi Biuxi Cu F 0 双曲型—波动方程
i 2 i 1
第一节 偏微分方程的有限差分解法
2.2 椭圆型偏微分方程的有限差分解法 用差分法将连续问题离散化的步骤： （1）首先对求解区域作网格剖分，用有限个网格节点代替连 续区域；
U xi 1, y j U xi , y j U xi 1, y j U xi , y j
U xi 1, y j U xi 1, y j 2U xi , y j
U xi 1 , y j 2U xi , y j U xi 1 , y j h12
U f x, y ,
hUij Ui 1, j 2Uij Ui 1, j h
2 1
x, y
Ui , j 1 2Uij Ui , j 1 h
2 2
fij
2 截断误差为 O h
第一节 偏微分方程的有限差分解法
1.2 椭圆型偏微分方程的有限差分解法 1.2.1 差分格式—— 九点差分格式
2 1 U xi , y j 2 U xi 1 , y j U xi , y j h1 h1 2 x 2! x 3 4 5 U x , y U x , y 1 1 1 U xi , y j 5 i j i j 3 4 6 h h h O h 1 1 1 1 3 4 5 3! x 4! x 5! x 2 1 U xi , y j 2 U xi 1 , y j U xi , y j h1 h1 2 x 2! x 3 4 5 1 U xi , y j 3 1 U xi , y j 4 1 U xi , y j 5 6 h h h O h 1 1 1 1 3! x3 4! x 4 5! x 5
地球物理计算方法
地球物理计算方法
课堂情况反馈
复习 上节课讲了些什么？
问题 （数值积分问题）
数值方法 （高斯积分公式）
复习
适当地选取求积节点（求系数Ak，待求节点xk) ，使求积公式具 有2n-1次代数精度（注意： n个求积节点）。
1
n
f (x)dx
1
Ak f (xk )
k 1
具有2n-1高精度的求积公式为高斯公式，该待求节点xk为高斯点。
效数字）
解：中点公式：
f 1 G h f 1 h f 1 h e1h e1h
2h
2h
h
G(h)
G1(h)
0.8
3.01765
G2(h)
G3(h)
0.4
2.79135
2.715917
0.2
2.73644
2.718137
2.718285
0.1
2.72281
2.719267
2.718276
2.71828
f
( x2 )
p2
(x0
th)
1 2
t
1t
2
f
(x0
)
t(t
2)
f
(x1)
1 2
t
t
1
f
( x2
)
对t求导，得
p
ห้องสมุดไป่ตู้
'2
(x0
th)
1 2h
2t
3
f
(x0
)
4(t
1)
f
(x1)
2t
1
f
( x2
)
t=0,1,2代入上式，得到三种三点公式
28
3、高阶导数公式 可以根据插值多项式构造：
课堂情况反馈
复习 上节课讲了些什么？
问题 （数值积分问题）
数值方法 （高斯积分公式）
复习
适当地选取求积节点（求系数Ak，待求节点xk) ，使求积公式具 有2n-1次代数精度（注意： n个求积节点）。
1
n
f (x)dx
1
Ak f (xk )
k 1
具有2n-1高精度的求积公式为高斯公式，该待求节点xk为高斯点。
效数字）
解：中点公式：
f 1 G h f 1 h f 1 h e1h e1h
2h
2h
h
G(h)
G1(h)
0.8
3.01765
G2(h)
G3(h)
0.4
2.79135
2.715917
0.2
2.73644
2.718137
2.718285
0.1
2.72281
2.719267
2.718276
2.71828
f
( x2 )
p2
(x0
th)
1 2
t
1t
2
f
(x0
)
t(t
2)
f
(x1)
1 2
t
t
1
f
( x2
)
对t求导，得
p
ห้องสมุดไป่ตู้
'2
(x0
th)
1 2h
2t
3
f
(x0
)
4(t
1)
f
(x1)
2t
1
f
( x2
)
t=0,1,2代入上式，得到三种三点公式
28
3、高阶导数公式 可以根据插值多项式构造：
第1章 地球物理中的计算问题概述
第五节 地球物理中的电磁场方程
电磁勘探从理论基础并结合异常的原因分为： （ 2 ）感应类电法，观测和利用的是由电磁感应作用产生的 异常，如瞬变电磁法、大地电磁测深法、频率测深法、无 线电波法等，它们满足麦克斯韦方程，且又是一种波场， 与地震波场在某些方法有相通之处。这种方法称为电磁法。
瞬变电磁法—利用不接地回线或接地线源向地下发射一次脉冲 磁场，在一次脉冲磁场间歇期间利用线圈或接地电极观测地下 介质中引起的二次感应涡流场，从而探测介质电阻率的一种方 法。 频率测深法—指频率在几十到几万Hz的音频范围内，通过改变 交变磁场频率的办法探测岩层电阻率随深度的变化以了解地质 构造和找矿的一种人工场源电磁法。 无线电波——频率在几十万赫至几十兆赫的电磁波。
理问题为出发点和归宿点。 从地球物理问题提出的方程、公式和求解方法与计算方
法都是比较独特的，不但与现成的数学方法不完全一致，而
且有时很难对它们的适定性、收敛性和稳定性做出预先的研 究。常常是先应用，从应用的效果来判断模型、方法和计算
过程的可行性。
##
第二节 地球物理中的引力位方程
P x, y, z
功与运动所经过的路径无关而只决定于运动的起始 点和终点的位置。这个性质叫力场的保守性。
第二节 地球物理中的引力位方程
对于保守力场可引入一个标量函数来描述。固定一点Q，对空 间任意点P（M0点除外）赋予一个标量V(P)—单位质点从Q点移 动到点P引力所作的功： Q
L
P x, y, z
Q处于无穷远处 F
许多以前无法计算的问题得以解决。这样才逐渐形成了计
算地球物理学的分支学科。 形成的开始时期大约在1960年代初期。
第一节 绪论
1.2 计算地球物理学的研究内容 正问题 物理模拟、解析解、数值计算 解析法——分离变数法、积分变换方法、 Green 函数方法、 变分法、对于二维和三维 Laplace 方程的边值问题，也还 可以将解表示为特殊的积分公式、保角变换，对于双曲型 方程的定解问题，也存在一些特殊的解法，例如平均值法， 降维法。
地球物理计算方法课件：第二章_数值积分 3
R
I
S n
ba 180
h 4
2
f
(4) ( )
b a h4 2880
f
(4) ( ),
(a,b).
n , h4 0, R 0
具有相应的收敛性和稳定性.
10
复化求积公式的截断误差
复化梯形公式的截断误差：
RT ( f )
b
f (x)dx T
a
n
h3 12
nf
(
)
ba 12
h2
f
(
)
I S2n 1 I Sn 16
移项整理，知
I
16 15
S
2n
1 15
Sn
16
1
这个值正好是科特斯公式的结果。也可写为 Cn 15 S2n 15 Sn
Cn (1 )S2n Sn
1
15
3、Romberg公式
重复同样的手续，依据柯特斯法的误差公式可进一步导出下列 龙贝格公式：
Rn
64 63 C2n
1 63 Cn
这个公式也可写为
Rn (1 )C2n Cn
1
63
18
一般递推公式：
k次二分得到得到的m次加速值，龙贝格一般公式
T (k) m
4m 4m
1
T (k m1
1)
1 4m 1
T( m1
k
)
.k
1,2,....
19
4、 Romberg算法的实现
计算流程见P70,实现步骤：
① 用梯形公式计算积分近似值
xk a kh(k 0,1,, n),
由定积分性质知
I
b
n1
f (x)dx
计算地球物理课件 第2章 地球物理中常用数值解法的基本原理-2
在满足边值条件的一切可能位置中，使位能取最小者。设弦处于
某一位置 u x ，计算它的位能(W 1 T 2, 为伸缩率):
2
应变位能 外力作功
第二节 偏微分方程的有限元解法
2.2 两点边值问题——弦的平衡 总位能
根据极小位能原理， u* u* x 是下列变分问题的解：
J
u*
min u
J
u
第二节 偏微分方程的有限元解法
几个概念
测度：有界开集和有界闭集的测度是区间长度的直接
推广。
E 是有界集 存在常数 M ，使对任意的 x (x1, x2, , xn ) E ，都有| xi | M (i 1, 2, , n) .
有界集
E
的外测度——
m*E
inf
Ii
,
Ii E ， inf 表
i 1
i 1
示最左边的意思。
有界集 E 的内测度——有界集 E 所包含的一切有界闭
第二节 偏微分方程的有限元解法
有限元法的基本问题可归纳为： （1）把问题转化成变分形式； （2）选定单元的形状，对求解域作剖分； （3）构造基函数或单元形状函数； （4）形成有限元方程（Ritz-Galerkin方程）； （5）提供有限元方程的有效解法； （6）收敛性及误差估计。
第二节 偏微分方程的有限元解法
x y z x y z （结合律）；
2）对任何 , k ， x, y X ，定义数乘，即 x X ，且满足
x x x ; x x ; x y x y ; 1 x x ;
3）在 X 中存在零元素，记为“0”，它满足
x0 x 4）对每个 x X ，存在 x 的加法逆元素，记为“-x” X ，使 x x 0
正定：设A是n阶实系 数对称矩阵，如果对 任何非零向量x都有 xTAx>0，就称A正定。
某一位置 u x ，计算它的位能(W 1 T 2, 为伸缩率):
2
应变位能 外力作功
第二节 偏微分方程的有限元解法
2.2 两点边值问题——弦的平衡 总位能
根据极小位能原理， u* u* x 是下列变分问题的解：
J
u*
min u
J
u
第二节 偏微分方程的有限元解法
几个概念
测度：有界开集和有界闭集的测度是区间长度的直接
推广。
E 是有界集 存在常数 M ，使对任意的 x (x1, x2, , xn ) E ，都有| xi | M (i 1, 2, , n) .
有界集
E
的外测度——
m*E
inf
Ii
,
Ii E ， inf 表
i 1
i 1
示最左边的意思。
有界集 E 的内测度——有界集 E 所包含的一切有界闭
第二节 偏微分方程的有限元解法
有限元法的基本问题可归纳为： （1）把问题转化成变分形式； （2）选定单元的形状，对求解域作剖分； （3）构造基函数或单元形状函数； （4）形成有限元方程（Ritz-Galerkin方程）； （5）提供有限元方程的有效解法； （6）收敛性及误差估计。
第二节 偏微分方程的有限元解法
x y z x y z （结合律）；
2）对任何 , k ， x, y X ，定义数乘，即 x X ，且满足
x x x ; x x ; x y x y ; 1 x x ;
3）在 X 中存在零元素，记为“0”，它满足
x0 x 4）对每个 x X ，存在 x 的加法逆元素，记为“-x” X ，使 x x 0
正定：设A是n阶实系 数对称矩阵，如果对 任何非零向量x都有 xTAx>0，就称A正定。
地球物理方法
地球物理方法地球物理方法:利用各种仪器在地表观测地壳上的各种物理现象，从而推断、了解地下的地质构造特点，寻找可能的储油构造。
它是一种间接找油的方法。
特点：精度和成本均高于地质法，但低于钻探方法。
地震勘探：就是利用人工方法激发的地震波(弹性波)，研究地震波在地层中传播的规律，以查明地下的地质构造，从而来确定矿藏（包括油气、矿石、水、地热资源等）等的位置，以及获得工程地质信息。
地震波运动学：研究在地震波传播过程中的地震波波前的空间位置与其传播时间的关系，即研究波的传播规律，以及这种时空关系与地下地质构造的关系。
地震波动力学：研究地震波在传播过程中波形、振幅、频率、相位等特征的及其变化规律，以及这些变化规律与地下的地层结构，岩石性质及流体性质之间存在的联系。
地震波：是一种在岩层中传播的，频率较低(与天然地震的频率相近)的波，弹性波在岩层中传播的一种通俗说法。
地震波由一个震源激发。
地震子波：爆炸产生的是一个延续时间很短的尖脉冲，这一尖脉冲造成破坏圈、塑性带，最后使离震源较远的介质产生弹性形变，形成地震波，地震波向外传播一定距离后，波形逐渐稳定，成为一个具有2-3个相位（极值）、延续时间60-100毫秒的地震波，称为地震子波。
地震子波看作组成一道地震记录的基本元素。
波前：振动刚开始与静止时的分界面，即刚要开始振动的那一时刻。
射线：是用来描述波的传播路线的一种表示。
在一定条件下，认为波及其能量是沿着一条“路径”从波源传到所观测的一点P。
这是一条假想的路径，也叫波线。
射线总是与波阵面垂直，波动经过每一点都可以设想有这么一条波线。
振动图和波剖面：某点振动随时间的变化的曲线称为振动曲线，也称振动图。
地震勘探中，沿测线画出的波形曲线，也称波剖面。
折射波：当入射波大于临界角时，出现滑行波和全反射。
在分界面上的滑行波有另一种特性，即会影响第一界面，并激发新的波。
在地震勘探中，由滑行波引起的波叫折射波，也叫做首波。
入射波以临界角或大于临界角入射高速介质所产生的波滑行波：由透射定律可知，如果V2>V1 ，即sinθ2 > sinθ1 ，θ2 > θ1。
地球物理计算方法课件：绪论2
(x)
1 10l1
2(1 1)
则x至少具有l位有效数字（即至少精确到它的第l位）。
例: 当用3.1416来表示 的近似值时，它的相对误差是多少？
例 为了使 x 20 的近似值x的相对误差小于0.1%，问至少
取几位有效数字？
误差的传播与估计
误差的传播
两个近似数x1与x2，其误差限分别为ε(x1) 及ε(x2)， 它们进行加、减、乘、除运算得到的误差限分别为
解为 x1=-6.222…, x2=38.25…, x3=-33.65…,
误差来源
用计算机解决科学计算问题(地球物理)的过程如下：
修修改改处处理理
地地质问质题 问题
物物理理模 模型型
数数学学模 模型型
数数值值计 算计算
否否
计足计要算算求结果结是否果满 是否满足
要求
是是 地地质质解 解释释
模型 误差
Cp越大，病态越严重.
避免误差危害的若干原则
1. 要避免除数绝对值远远小于被除数绝对值的除法 用绝对值小的数作除数舍入误差会增大，如计算 x
y
若0 y ， 则x 可能对计算结果带来严重影响，应 尽量避免.
2. 要避免两相近数相减
在数值计算中两相近数相减有效数字会严重损失. 例如， x 532.65, 都y 具5有32五.5位2有效数字， 但 x y 只0.有13两位有效数字. 这说明必须尽量避免出现这类运算. 最好是改变计算方法， 防止这种现象产生.
a1 为0整, m数，取l
x x * 1 10ml. 2
0ln
有效数字与绝对误差限的关系
对 取3位和5位有效数字，他们的误差限是多少？
3* 0.002, 5* 0.000008,
地球物理方法技术-电磁法
2.3 电磁法(EM)
2.3.4 电磁法的应用实例
在该地主要是寻找裂隙水，断裂探测为主要目标。在该区布置了勘探剖面，EH4电磁成像电阻 率断面如图2，发现了两个有利的断裂，在F1断裂处实施钻探，在580m处获得高产水。
2.3 电磁法(EM)
2.3.4 电磁法的应用实例
上团
21 21 21 9 5 1
2.3 电磁法(EM)
2.3.3 电磁法的异常特征 3) 甚低频法
f: 15500—22300Hz
2.3 电磁法(EM)
2.3.3 电磁法的异常特征 4) 音频大地电流法
f: 1—20KHz
2.3 电磁法(EM)
2.3.4 电磁法的应用实例 1)联合使用Slingram和常规VLF法勘探地下水
垂直线框
-----
2.3 电磁法(EM)
2.3.3 电磁法的异常特征 2 电 磁 偶 极 源 虚 分 量 振 幅 法 ) )
水平线框
-----
2.3 电磁法(EM)
2.3.3 电磁法的异常特征 2 电 磁 偶 极 源 虚 分 量 振 幅 法 ) )
旁线(XZ)
-----
2.3 电磁法(EM)
2.3.3 电磁法的异常特征 2 电 磁 偶 极 源 ) ) 水 平 线 圈 法
Specialist Training Course on EMD
地球物理方法技术
中国地质大学(北京)地球物理学院 杨 进
地球物理方法技术
第1章. 第2章. 第3章. 第4章. 第5章.
地球物理方法绪论 传统的地球物理方法 地球物理新方法、新技术 地质问题中地球物理方法的策略 地球物理勘探实例
地球物理方法技术 第2章. 传统的地球物理方法 2.1 电阻率法(DC) 2.2 激发极化法(IP) 2.3 电磁法 (EM) 2.4 大地电磁法 (MT) 2.5 磁法和重力(MM,GM) 2.6 浅层地震反射法 2.7 测井
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
算法流程图
计算方法及地球物理数值计算 算法设计 误差分析 Matlab计算方法基础
误差分析必要性
例：求解方程 准确解为:
x2 (105 1)x 105 0
x1 105, x2 1;
数值近似解（5位有效数字）：
x1 105, x2 0;
避免大数”吃掉”小数 ;
例：方程组的解
x1
第一步：计算f (x)在端点处的值f (a), f (b).
第二步：计算f
( x)在中点处的值x1
a
b 2
,
f
(a
2
b)
第三步：若f (x ) 0
则取 x* x1 a
b 2
若f
(x
)
f
(a
)
0,则取 a1
a, b1
ab 2
若f
(x
)fBiblioteka (b ）0,则取 a1
a
2
b
，b1
b
重复步骤2和步骤3.
1 2
x2
1 3
x3
11 6
1
2
x1
1 3
x2
1 4
x3
13 12
1 3
x1
1 4
x2
1 5
x3
47 60
准确解为 x1=x2=x3=1
系数舍入近似为
x1 0.50x2 0.33x3 1.8 0.50x1 0.33x2 0.25x3 1.1 0.33x1 0.25x2 0.20x3 0.78
解为 x1=-6.222…, x2=38.25…, x3=-33.65…,
误差来源
用计算机解决科学计算问题(地球物理)的过程如下：
修修改改处处理理
地地质问质题 问题
物物理理模 模型型
数数学学模 模型型
数数值值计 算计算
否否
计足计要算算求结果结是否果满 是否满足
要求
是是 地地质质解 解释释
模型 误差
e x x*
为近似值的绝对误差，简称误差.
误差 e可正可负， 当绝对误差为正时近似值偏大，叫强近似值 当绝对误差为负时近似值偏小，叫弱近似值
地球物理计算方法
张致付 地球物理与信息技术学院
绪论
复习
计算方法及地球物理数值计算 算法设计 误差分析 Matlab计算方法基础
计算方法
计算方法
公式、算法（方法） 理论分析（收敛性, 稳定性,
误差分析等）
地球物理数值计算
•数据处理
（插值或拟合、FFT、矩阵求解、数值积分）
•数值模拟正演
（微分方程求解）
写成递推公式:
vv0k
vk 1x ank an
(k 1, 2, , n)
秦九韶算法,只做n次乘法和n次 加法,程序实现简单
算法流程图
3、迭代校正-逼近
设有非线性方程
f (x) 0
其中，f(x)在[a,b]上连续函数且f(a)f(b)<0，假设方程于[a,b] 内仅有一个实根。
求非线性方程f(x)=0的根，二分法分为两步： ➢确定方程的有根区间 ➢计算根的近似值（二分法）
数值计算是用计算机求解， 计算机不能做无限运算、 不能 准确表达一个实数。受计算机字长限制产生舍入误差，少量的 舍入误差一般是微不足道的， 而计算机上完成成千上万次运算 ，舍入误差的积累就大了，可能导致错误的结果。
例如：用3.14159代替
绝对误差和相对误差
误差：设x*为准确值， x为 x*的一个近似值，称
观测误差: 由物理量观测产生的误差，在数学建模和具
体计算过程中所用的数据往往是通过观察和测量得到的， 由于测量仪器精度的限制，这些数据一般是近似的。
物理量数据：温度、长度、电位、电流、速度等。
3、截断误差
截断误差：计算机只能完成有限次算术
运算和逻辑运算；是连续型数学问题向 离散型数值算法转化过程中产生的，用 有限过程近似无限过程的误差。这对无 穷过程进行截断所带来的误差。
算法设计必要性
例：
N阶线性方程组，求解
Crame 法则 Gauss列主元消去法
Crame 法则
需 要 多 少 次 乘 法 ?
Crame 法则
①一个 n 阶行列式有n!个乘积，每个乘积要作（n-1）次乘法 ②要计算n+1 个行列式 ③要 n 次除法
总计算量 n!(n-1)(n+1)+n，n=20 时要作9.7×1021次乘法。
截断 误差
舍入 误差
1、模型误差
➢ 首先要建立数学模型，
它是对被描述的实际问题进行抽象、简化而 得到的，因而是近似的。
模型误差：在建立数学模型过程中，要将复杂 的现象抽象归结为数学模型，建数学模型时加 上了一些限制，往往忽略了一些次要因素的影 响，与实际问题有误差。
实际问题 数学模型
2、观测误差
f (n1) ( )
(n 1)!
x n 1
在0与x之间。
如常用的函数Taylor展开：
ex =1 x x2 x3 ... 2! 3!
sin(x)=x x3 x5 x7 ... 3! 5! 7!
ln(1 x)=x x2 x3 x4 ... 2! 3! 4!
4、舍入误差
舍入误差：计算过程取有限位数
Gauss 列主元消去法
乘除法：n3/3+n2-n/3 只需做3060 次乘、除运算。
算法设计 1、科学计算可实现
例:证明二次方程
x2 2bx c 0
至多有两个不同的实根。
（1）反证法 （2）图解法 （3）公式法——算法
计算机只能做加、乘基 本算法
算法流程图
2、计算量的减少
例：多项式求和
设要求对给定x的求下列多项式的值
p x a0 a1x a2x2 an xn
直接求法
若直接计算akxk,再逐项相加，一共要做
n+(n-1)+…+2+1=n(n+1)/2次乘法和n 次加法。
秦九韶算法
如果将前n项提出x，则有
p(x) （an xn1 an1xn2 a1）x a0 ((an xn2 an1xn3 a2 )x a1）x a0 ((an x an1)x a2 )x a1)x a0
•物理参数反演计算
（线性方程组求解）
复习
计算方法及地球物理数值计算 算法设计 误差分析 Matlab计算方法基础
• 计算方法：研究如何通过计算机所能执行的基 本运算，求出各类实际问题的数值解。
• 算法：给出已知的量，通过给定的运算次序， 经过有限次的基本运算，得出所求的未知量的 解，这种完整的运算步骤称为算法。
数学模型 数值计算方法
截断误差（方法误差）:近似解与精确解之间的误差
例：用Taylor公式进行展开
Pn (x)
=f(0)+
f '(0) x 1!
f
''(0) x2 ... 2!
f (n) (0) xn n!
近似代替函数f(x),则数值方法的截断误差是：
Rn (x)
=f (x)-Pn (x)=