(二)河道砂体正演分析

浅层新近系地层以砂泥岩薄互层为特点。钻井统计表明: 80% 以上的储层单层厚度小于 10 m，即大部分砂层都属于薄层的范畴。由于地震分辨率的限制，无法分辨单一薄层，也就是说常规地震剖面上的一个同相轴是由若干个薄层反射叠加而成。浅层地层具有埋藏浅、压实程度低的特点，由于砂岩与泥岩速度往往相差较小，在有些地区馆陶组砂岩速度接近、甚至低于泥岩速度，这种状况造成了不同地区馆陶组砂岩地震反射特征的不同。因此，对于馆陶组的储层预测与描述工作而言，首先必须进行模型正演，借此来分析研究区砂岩储层在地震上的反射特征。建立起地质与地震之间的桥梁，为下一步定量储层预测打下坚实基础。

通过正演模拟可解决几个问题: ①识别有关储层的地震响应信息; ②在常规地震剖面上识别地层或岩性圈闭的可能性; ③建立属性参数与地层岩性参数的对应关系，为地震资料解释提供依据; ④利用地震属性参数确定地层、岩性参数，实现储层的定量化描述。

1地震正演模拟原理及实现方法

地震正演模拟以垂直入射线反射理论为算法基础。射线追踪的基本原理遵循惠更斯原理和费马原理。惠更斯原理的解释式是时间场特征方程式:

济阳坳陷北部馆陶组油气地质与勘探技术

式中: v2x，y，z为某一类型的波 (纵波或横波)在介质中给定的速度分布。

费马原理是地震波沿射线传播的时间和沿其他可能的途径传播所需的时间比较起来为最小。费马原理的数学表达是解下列积分极值的变分问题 (姜秀清，2002)。

济阳坳陷北部馆陶组油气地质与勘探技术

式中: s (r)=1/c (r)，为慢度; c (r)为介质的速度; μ 为射线参数; σ 为弧长，且有下式:

济阳坳陷北部馆陶组油气地质与勘探技术

对 (8-2)式变分求极值、化简可得二阶常微分方程，即

济阳坳陷北部馆陶组油气地质与勘探技术

式 (8-3)就是射线应满足的方程。

这两个原理的直接结果是射线追踪的斯奈尔定律:

济阳坳陷北部馆陶组油气地质与勘探技术

式中: p 为射线参量。

本次采用运动学 (射线)理论计算波的旅行时间，得出各种时距曲线; 考虑波的动力学特点，把界面的反射系数、透射系数、各种多次波、吸收损失、波前扩散等影响因素考虑进去，在计算出的地震响应中不只反映波的传播时间的特点，还可反映波的振幅等特点。

2正演模型设计

综合已知资料，在充分了解地下地质结构的基础上，建立速度-深度 (地层)模型，然后模拟实际资料所采用的野外观测方法，进行采集、处理，最后得到合成地震响应。将合成地震数据与实际地震数据进行比较，若偏差较大，则修正速度-深度模型，直到较为吻合为止，这就是地震模型的正演过程 (图 8-1)。

图 8-1 模型正演流程图

正演所得到的结果，一方面用于指导地震解释，另一方面速度、深度参数可用于油藏描述。对于薄层、岩性组合模型的正演，关键是分析地震响应的波形变化，应与模型技术的反演过程结合起来，从实际地震道中提取子波，用于正演模拟 (姜秀清，2002)。

模型的正演及其修正过程实际上是地震资料-地质解释反复认识和深化的过程，其目的是寻求合理地质解释，减少地震解释的多解性。

3观测系统

建立模型后，就要设计观测系统，运用射线追踪来模拟地下地层。通过比较不同炮点(或接收点)的结构和地质组合模型，决定一个特殊的炮点 (或接收点)排列是否可以对地下目标成像。通过模型的射线道可以看到射线追踪地下每一层的过程。用任意射线追踪方法拾取任意炮点的位置，观察到可能产生的射线追踪矢量，反复修改排列参数，达到希望追踪的结果，设计出最佳的炮点 (或接收点)的排列，为下一步地震处理提供良好的数据来源。

4地震处理

用模型软件变换和测试各种地质体，用岩性和流体含量来比较它们如何影响最终结果。需要进行以下工作:

(1)改变砂泥岩密度、速度;

(2)变换倾斜后层位深度;

(3)改变储层的流体特性，比较模拟结果。

利用实钻井统计的基础资料，根据本区地层结构和沉积微相特征，建立不同的模型系列，利用正演模拟技术对预测地质模型进行地震模拟，从而分析储层的地震响应特征并与实际地震剖面对比，不断地调整地质模型，使模拟地震剖面与实际地震剖面最大的相似，因而能够大体掌握井点以外目的层的展布范围，在此认识的基础上，再进行地震反演、地震属性分析，使预测结果更加真实、可信。

5油藏模型响应特征分析

(1)楔状体模型

根据该区实际统计的砂泥岩速度设计的楔状体模型表明，振幅、频率随薄层厚度变化具有一定的变化规律 (图 8-2)，储层厚度大于调谐厚度 λ/4 (155 m)时，储层的顶底反射随厚度增加逐渐分开，形成两个反射同相轴，且反射振幅值随厚度的增加逐渐趋于稳定，反射波峰与波谷时间差反映了砂体的真实厚度，称为时间可分辨区; 若改变砂泥岩的速度差，可以看到只是砂泥岩的反射界面强度改变，而波形开始畸变的点不变。当砂体厚度小于 λ/4 时，砂体的顶底反射叠加在一起，其波峰与波谷的时间差基本保持不变，趋近于一个常数，时差值已不能反映储层厚度的信息，这时砂体反射的相对振幅值与厚度呈近似线性关系，即地震正演记录的振幅随储层厚度的减小而明显减小，在没有砂体的地方地震正演记录出现 “相位上拉”现象。因此可由振幅预测厚度，称为时间不可分辨的振幅可分辨区，当砂体厚度等于 λ/4 时，砂体顶底反射相干加强，振幅最大，λ/4 称为调谐厚度。因此，当储层厚度大于调谐厚度时，不能直接描述单个同相轴来解释砂体，可以直接根据地震波的反射旅行时计算储层厚度; 当储层厚度小于调谐厚度时，则出现一个同相轴，可以直接利用该同相轴的波峰或波谷描述单砂体，而且可用振幅能量按其线性关系计算储层的厚度。

图 8-2 不同砂泥岩条件下楔状体模型及地震响应特征

对总厚度小于 λ/4 的地质模型及其地震响应特征分析表明，由厚度不同的砂岩薄层组成的砂泥岩薄互层虽然内部排列方式不同，但其地震响应特征大致相同。

(2)不同沉积微相地震反射特征

济阳坳陷北部馆陶组沉积砂体包括河道型、砂坝型、漫滩型等储集体类型、 (高速、低速)两种围岩介质。由于砂质泥岩和泥质砂岩速度高于砂岩，因此地震反射较为丰富，且不同区带存在一定的差异。

1)正演模型设计。由于曲流河沉积环境中砂岩百分含量小，为泥包砂的特征，泥岩隔层厚; 通过对该区河道发育区的地震资料统计，频率在 33Hz 左右。在模型设计过程中，储层之间相对不会出现干扰 (即单储层模型设计特征可反映网状河的特征)，因此我们对单套储层进行研究。在模型设计过程中主要讨论曲流河不同沉积微相的沉积特征。不同沉积微相的储层及非储层具有不同的岩石物理特征，在曲流河沉积时期受各种因素的影响，在纵向上和横向上形成不同的组合与接触关系，由此产生不同的反射特征。

测井相就是与沉积相及储集层特征相联系的测井响应特征的总和，通过对测井曲线进行研究可确定不同的地质体类型。目前应用最广的是自然电位曲线和视电阻率曲线，特别是自然电位曲线，自然电位曲线能较好的反映垂向上的沉积演化。在模型设计过程中，主要根据自然电位曲线特征设计不同储盖层的组合方式。

自然电位是地层扩散-吸附电势、过滤电势以及氧化还原电势的综合响应。扩散-吸附电势取决于地层水和钻井液之间离子的浓度差、岩层中泥岩含量，以及受粒度和分选控制的孔吼半径大小。地层渗透性决定过滤电势，渗透性又与粒度、分选和泥质含量有关。砂、泥岩层的氧化还原电势取决于水动力条件的强弱。因此认为，由三种电势构成的自然电位主要受粒度、分选和泥质含量的控制，而它们又受沉积时水动力能量和物源供应条件的制约，所以自然电位曲线的变化能反映沉积环境的变迁。当砂岩自下而上粒度变小、分选变差、泥质含量增多时，则意味着水动力能量变弱，物源供应减少，因而自然电位曲线幅度也向上变小。自然电位曲线的幅度、形态、顶底接触关系、曲线光滑程度代表沉积过程的沉积响应。根据变化可分析出不同的组合关系。

模型设计过程中，主要根据几种常见的自然电位曲线 (钟形、漏斗形、箱形)结合其他特征进行设计。

分析认为河流相砂体沉积时期主要发育两种围岩: 一种是速度较高的围岩 (即砂质泥岩)，泥岩相对不发育; 另一种是低速围岩，泥岩发育，而砂质泥岩不发育。通过对不同地区的速度统计，具有泥岩、砂岩、泥质砂岩速度依次增大的结论。砂岩速度一般为2500 m/s，泥岩速度一般为2200m/s，泥质砂岩与砂质泥岩一般为2600m/s; 含油砂岩一般为2450m/s，含气砂岩一般为1800m/s，岩石中流体为含油水层或油水同层对岩石速度影响不明显。密度统计结果为: 泥岩密度20 ～21g/cm3，砂岩密度 195 ～ 215 g/cm3，泥质砂岩密度 205 ～233 g / cm3，含气砂岩密度 185 ～ 195 g/cm3，含油砂岩密度 195 ～ 203 g/cm3。从小到大依次顺序为: 含气砂岩、含油砂岩、砂岩、泥岩、泥质砂岩。该统计结果是我们建立模型的基础。

图8-3 钟形相序沉积体正演模型的建立

2)正演特征分析

A钟形: 常出现于河道型储集体中。反映在该沉积时期水流能量逐渐减弱，物源供应不断减少，在物源丰富和水动力作用强的条件下，被充分淘洗后的均质沉积，底部加速式渐变，在冲刷面下部有原先滞留的沉积砂; 顶部相对匀速渐变，代表匀速的能量减退;在该沉积时期，物源逐渐减少，水流冲刷、簸扬能力逐渐减弱，砂体改造彻底，分选较好，为一套较好的储集层。模型设计中，储层底部为分选较好的砂子，由于水流能量减弱和物源供应的不断减少，沉积后期为一渐变的分选逐渐变差的储层，直到物性差的围岩或极少发育的泥岩覆盖其上。

a顶部覆盖高速围岩，底部接触高速围岩: 为单轨强反射，两边弱上拉。储层底面为一正极性反射特征反射，略大于反射的最大振幅处 (图 8-3a)。

b顶部覆盖低速泥岩，底部接触高速围岩: 为双轨中强反射。分别在顶部低速泥岩与变速围岩的接触面和底部高速围岩与储层底部的接触面产生一套正极性反射特征的反射，低速泥岩与变速围岩的接触面对应产生反射的最大振幅处，底部高速围岩与储层底部的接触面对应产生反射并不对应最大振幅处 (图 8-3b)。

c底部接触低速泥岩，顶部覆盖高速围岩: 储层顶面为弱反射或无反射。只在低速泥岩与围岩产生一正极性反射界面，储层与其他接触面不产生反射界面 (图 8-3c)。

d底部接触低速泥岩，顶部覆盖低速泥岩: 为强反射。对应于顶部低速泥岩与变速围岩的接触面，且该接触面对应于产生反射的最大振幅处 (图 8-3d)。

B漏斗形: 常出现于砂坝型储集体中。与钟形的形态相反，反映在该沉积时期向上水流能量加强，物源供应不断增加，分选逐步变好。砂体在沉积初期物源供应不足，沉积后期为一渐变的分选逐渐变好的储层，直到物性差的围岩或极少发育的泥岩覆盖其上 (图 8-4)。

模型设计中，分两种情况进行讨论。储层顶部底部均接触高速围岩: 正演结果为一套强正极性反射，在储层顶面与高速围岩的接触面上对应产生一负极性反射且在反射的最大振幅处 (图 8-4a)。

在接触的围岩顶部覆盖一层薄低速泥岩: 产生两套正极性反射特征反射且下部反射强于上部反射 (图 8-4b)。

C箱形: 沉积过程中物源丰富，水流能量加强且水动力条件稳定，沉积初期存在冲刷面，沉积后期物源供应的突然中断，水动力突然变小，可能是截弯取直的废弃河道。

图 8-4 漏斗形相序沉积体正演模型的建立

箱形模型比较简单，其在储层顶底均接触突变的高速围岩: 正演结果储层的顶部为一负极性反射特征反射且在反射的最大振幅处，底部出现一正反射特征的反射，但不在反射的最大振幅处 (图 8-5a)。

箱型模型，其在储层顶低均接触突变的低速围岩: 正演结果储层的顶部为一正极性反射特征反射，且砂体顶面在反射的最大振幅处 (图 8-5b)。

从设计的几个模型所正演的结果可以看到，与储层接触的围岩是产生不同反射的主要因素，尤其是低速泥岩的出现直接影响到反射特征的改变。我们在地震剖面上看到的反射并不一定是储层的顶面或是底面反射。

3)薄互层沉积砂体正演模型地震反射特征

由于辫状河、曲流河、三角洲、湖泊沉积等沉积类型，以及沉积微相类型在纵向上的快速变化及交互叠置，砂泥薄互层也是济阳坳陷北部馆陶组常见的储盖组合。在频率33Hz 左右的情况下，模型设计过程中储层之间可能会由于隔层太薄产生干扰，因此我们对两套储层进行研究，讨论隔层的厚度对储层反射产生的影响，同时对不同砂泥岩比导致的不同的反射特征进行研究。

图 8-5 箱形正演模型的构成及其特征

模型 1: 在相同环境相同深度下分别设计两套厚度相同的储层，避免由于环境和深度对砂体反射产生的影响，单纯考虑隔层厚度因素。模型设计中，两砂体之间分别间隔 5 m、10 m、15 m、20 m、25 m 不同厚度的围岩隔层，从正演结果可以看到，当隔层为 5 m 时，在第一套储层处出现一套正极性反射; 隔层为 10 m 时，也在第一套储层处出现一套正极性反射，但反射强度明显弱于隔层为 5 m 时的反射; 当隔层大于10 m 时，在第一、二套储层处均出现一套正极性反射，随隔层厚度增加，反射强度增大 (图 8-6)。

模型 2: 相同砂泥岩比，不同的储盖组合产生不同的地震响应特征 (图 8-7)。

模型 3: 低砂泥比的砂泥岩薄互层，横向上砂体叠置连片，但砂体之间往往有隔层分开。针对这两类不同沉积特征的曲流河砂体分别设计正演模型 (图 8-8)。

图 8-8 为具有纵向叠置河道型砂体的地震响应波阻抗模型正演记录。可以看出，在围岩速度大于砂岩速度时，地震正演记录的振幅随储层厚度的减小而明显减弱，在河道型砂体叠置的地方地震正演记录出现 “振幅增强”现象。所以，一般情况下，在没有其他条件证明是叠置河道砂体时，解释人员自然也将其 “综合”为一条 (大的)河道砂体进行解释了。

图 8-6 不同厚度围岩间隔正演模型图

图 8-7 不同储盖组合正演模型图

图 8-8 砂体叠置河道砂体正演模型

D多井复合模型

图8-9、8-10 是根据实际钻井的岩性、地层及速度特征设计的地质模型及地震响应特征。

地震正演模拟结果表明，馆陶组河道型砂岩大部分具有中强振幅反射特征，只有上覆围岩阻抗值低于砂岩阻抗值时，砂体表现为弱反射。砂体集中在700 ～1400ms 之间，河道内反射杂乱，不连续; 外围呈较平行、连续反射。馆陶组 1 砂组: 呈强短轴反射，横向变化快。河道砂体表现为强反射特征，地震响应模式为中强振幅、低频、低阻抗，延伸范围较短，有的有下拉现象、有的有上凸现象，与围岩有明显的突变接触关系。这种强反射是一种组合反射效应。因为河道砂体有效厚度平均仅 5 m，尤其是小于5m 的砂岩厚度在常规地震资料中是无法产生强反射波的 (杜劲松，2004)。强反射的形成是河道砂岩与上覆高速的砂质泥岩或灰质泥岩组合的反射结果。馆陶组 2 砂组: 反射轴连续。砂体的平面展布为典型曲流河形态，有的主要为平面上形态不规则的点砂坝沉积，中强振幅，但幅值变化较大。

(3)河道型砂体分辨率问题

砂体间距小于一个采样间隔的砂体，在主频较高时其反射仍然连续，这容易造成解释的误区，造成钻井部署的失误; 而近距离叠置的砂体，低频时无法区分，主频较高时可以区分开。用小道距、高分辨率采集资料可以避免这种现象 (图 8-11)。

图 8-9 济阳坳陷北部馆陶组河道型砂体地质模型

图 8-10 济阳坳陷北部馆陶组河道型砂体正演模拟剖面

图 8-11 砂体分辨率正演模型图 (从上往下编号依次为 a、b、c)

不同频率下，地震响应特征有差异; 不同厚度砂体的地震优势频率和分辨率是不同的，既有强反射，又有弱反射 (图 8-12a、b、c)。

图 8-12 砂体模型及合成记录

由图 8-13 可见，当信噪比 (S/N)小于 3 时，薄层叠合砂体反射难以识别，当信噪比小于 2 时，厚层砂体反射边界模糊不清。可见，信噪比影响砂体描述。

图8－13 信噪比对地震分辨率的影响

摘 要：介绍了电磁学计算方法的研究进展和状态，对几种富有代表性的算法做了介绍，并比较了各自的优势和不足，包括矩量法、有限元法、时域有限差分方法以及复射线方法等。

关键词：矩量法；有限元法；时域有限差分方法；复射线方法

1 引 言

1864年Maxwell在前人的理论（高斯定律、安培定律、法拉第定律和自由磁极不存在）和实验的基础上建立了统一的电磁场理论，并用数学模型揭示了自然界一切宏观电磁现象所遵循的普遍规律，这就是著名的Maxwell方程。在11种可分离变量坐标系求解Maxwell方程组或者其退化形式，最后得到解析解。这种方法可以得到问题的准确解，而且效率也比较高，但是适用范围太窄，只能求解具有规则边界的简单问题。对于不规则形状或者任意形状边界则需要比较高的数学技巧，甚至无法求得解析解。20世纪60年代以来，随着电子计算机技术的发展，一些电磁场的数值计算方法发展起来，并得到广泛地应用，相对于经典电磁理论而言，数值方法受边界形状的约束大为减少，可以解决各种类型的复杂问题。但各种数值计算方法都有优缺点，一个复杂的问题往往难以依靠一种单一方法解决，常需要将多种方法结合起来，互相取长补短，因此混和方法日益受到人们的重视。

本文综述了国内外计算电磁学的发展状况，对常用的电磁计算方法做了分类。

2 电磁场数值方法的分类

电磁学问题的数值求解方法可分为时域和频域2大类。频域技术主要有矩量法、有限差分方法等，频域技术发展得比较早，也比较成熟。时域法主要有时域差分技术。时域法的引入是基于计算效率的考虑，某些问题在时域中讨论起来计算量要小。例如求解目标对冲激脉冲的早期响应时，频域法必须在很大的带宽内进行多次采样计算，然后做傅里叶反变换才能求得解答，计算精度受到采样点的影响。若有非线性部分随时间变化，采用时域法更加直接。另外还有一些高频方法，如GTD，UTD和射线理论。

从求解方程的形式看，可以分为积分方程法（IE）和微分方程法（DE）。IE和DE相比，有如下特点：IE法的求解区域维数比DE法少一维，误差限于求解区域的边界，故精度高；IE法适合求无限域问题，DE法此时会遇到网格截断问题；IE法产生的矩阵是满的，阶数小，DE法所产生的是稀疏矩阵，但阶数大；IE法难以处理非均匀、非线性和时变媒质问题，DE法可直接用于这类问题〔1〕。

3 几种典型方法的介绍

有限元方法是在20世纪40年代被提出，在50年代用于飞机设计。后来这种方法得到发展并被非常广泛地应用于结构分析问题中。目前，作为广泛应用于工程和数学问题的一种通用方法，有限元法已非常著名。

有限元法是以变分原理为基础的一种数值计算方法。其定解问题为：

应用变分原理，把所要求解的边值问题转化为相应的变分问题，利用对区域D的剖分、插值，离散化变分问题为普通多元函数的极值问题，进而得到一组多元的代数方程组，求解代数方程组就可以得到所求边值问题的数值解。一般要经过如下步骤：

①给出与待求边值问题相应的泛函及其变分问题。

②剖分场域D，并选出相应的插值函数。

③将变分问题离散化为一种多元函数的极值问题，得到如下一组代数方程组：

其中：Kij为系数（刚度）矩阵；Xi为离散点的插值。

④选择合适的代数解法解式（2），即可得到待求边值问题的数值解Xi（i＝1，2，…，N）

（2）矩量法

很多电磁场问题的分析都归结为这样一个算子方程〔2〕：

L（f）＝g（3）其中：L是线性算子，f是未知的场或其他响应，g是已知的源或激励。

在通常的情况下，这个方程是矢量方程（二维或三维的）。如果f能有方程解出，则是一个精确的解析解，大多数情况下，不能得到f的解析形式，只能通过数值方法进行预估。令f在L的定义域内被展开为某基函数系f1，f2，f3，…，fn的线性组合：

其中：an是展开系数，fn为展开函数或基函数。

对于精确解式（2）通畅是无限项之和，且形成一个基函数的完备集，对近似解，将式 （2）带入式（1），再应用算子L的线性，便可以得到：

m＝1，2，3，…

此方程组可写成矩阵形式f，以解出f。矩量法就是这样一种将算子方程转化为矩阵方程的一种离散方法。

在电磁散射问题中，散射体的特征尺度与波长之比是一个很重要的参数。他决定了具体应用矩量法的途径。如果目标特征尺度可以与波长比较，则可以采用一般的矩量法；如果目标很大而特征尺度又包括了一个很大的范围，那么就需要选择一个合适的离散方式和离散基函数。受计算机内存和计算速度影响，有些二维和三维问题用矩量法求解是非常困难的，因为计算的存储量通常与N2或者N3成正比（N为离散点数），而且离散后出现病态矩阵也是一个难以解决的问题。这时需要较高的数学技巧，如采用小波展开，选取合适的小波基函数来降维等〔3〕。

（3）时域有限差分方法

时域有限差分（FDTD）是电磁场的一种时域计算方法。传统上电磁场的计算主要是在频域上进行的，这些年以来，时域计算方法也越来越受到重视。他已在很多方面显示出独特的优越性，尤其是在解决有关非均匀介质、任意形状和复杂结构的散射体以及辐射系统的电磁问题中更加突出。FDTD法直接求解依赖时间变量的麦克斯韦旋度方程，利用二阶精度的中心差分近似把旋度方程中的微分算符直接转换为差分形式，这样达到在一定体积内和一段时间上对连续电磁场的数据取样压缩。电场和磁场分量在空间被交叉放置，这样保证在介质边界处切向场分量的连续条件自然得到满足。在笛卡儿坐标系电场和磁场分量在网格单元中的位置是每一磁场分量由4个电场分量包围着，反之亦然。

这种电磁场的空间放置方法符合法拉第定律和安培定律的自然几何结构。因此FDTD算法是计算机在数据存储空间中对连续的实际电磁波的传播过程在时间进程上进行数字模拟。而在每一个网格点上各场分量的新值均仅依赖于该点在同一时间步的值及在该点周围邻近点其他场前半个时间步的值。这正是电磁场的感应原理。这些关系构成FDTD法的基本算式，通过逐个时间步对模拟区域各网格点的计算，在执行到适当的时间步数后，即可获得所需要的结果。

在上述算法中，时间增量Δt和空间增量Δx，Δy和Δz不是相互独立的，他们的取值必须满足一定的关系，以避免数值不稳定。这种不稳定表现为在解显式 差分方程时随着时间步的继续计算结果也将无限制的67增加。为了保证数值稳定性必须满足数值稳定条件：

其中：（对非均匀区域，应选c的最大值）〔4〕。

用差分方法对麦克斯韦方程的数值计算还会在网格中引起所模拟波模的色散，即在FDTD网格中数字波模的传播速度将随波长、在网格中的传播方向以及离散化的情况而改变。这种色散将导致非物理原因引起的脉冲波形的畸变、人为的各向异性及虚拟的绕射等，因此必须考虑数值色散问题。如果在模拟空间中采用大小不同的网格或包含不同的介质区域，这时网格尺寸与波长之比将是位置的函数，在不同网格或介质的交界面处将出现非物理的绕射和反射现象，对此也应该进行定量的研究，以保证正确估计FDTD算法的精度。在开放问题中电磁场将占据无限大空间，而由于计算机内存总是有限的，只能模拟有限空间，因此差分网格在某处必将截断，这就要求在网格截断处不引起波的明显反射，使对外传播的波就像在无限大空间中传播一样。这就是在截断处设置吸收边界条件，使传播到截断处的波被边界吸收而不产生反射，当然不可能达到完全没有反射，目前已创立的一些吸收边界条件可达到精度上的要求，如Mur所导出的吸收边界条件。

（4）复射线方法

复射线是用于求解波场传播和散射问题的一种高频近似方法。他根据几何光学理论和几何绕射理论的分析方法和计算公式，在解析延拓的复空间中求解复射线轨迹和场的振幅和相位，从而直接得出局部不均匀波（凋落波）的传播和散射规律〔5〕。复射线方法是包括复射线追踪、复射线近轴近似、复射线展开以及复绕射线等处理技术在内的一系列处理方法的统称。其共同特点在于：通过将射线参考点坐标延拓到复空间而建立了一个简单而统一的实空间中波束／射线束（Bundle ofrays）分析模型；通过费马原理及其延拓，由基于复射线追踪或复射线近轴近似的处理技术，构造了射线光学架构下有效的鞍点场描述方法等。例如，复射线追踪法将射线光学中使用的射线追踪方法和场强计算公式直接地解析延拓到复空间，利用延拓后的复费马原理进行复射线搜索，从而求出复射线轨迹和复射线场。这一方法的特点在于可以基于射线光学方法有效地描述空间中波束的传播，因此，提供了一类分析波束传播的简便方法。其不足之处是对每一个给定的观察点必须进行一次二维或四维的复射线轨迹搜索，这是一个十分花费时间的计算机迭代过程。

4 几种方法的比较和进展

将有限元法移植到电磁工程领域还是二十世纪六七十年代的事情，他比较新颖。有限元法的优点是适用于具有复杂边界形状或边界条件、含有复杂媒质的定解问题。这种方法的各个环节可以实现标准化，得到通用的计算程序，而且有较高的计算精度。但是这种方法的计算程序复杂冗长，由于他是区域性解法，分割的元素数和节点数较多，导致需要的初始数据复杂繁多，最终得到的方程组的元数很大，这使得计算时间长，而且对计算机本身的存储也提出了要求。对电磁学中的许多问题，有限元产生的是带状（如果适当地给节点编号的话）、稀疏阵（许多矩阵元素是0）。但是单独采用有限元法只能解决开域问题。用有限元法进行数值分析的第一步是对目标的离散，多年来人们一直在研究这个问题，试图找到一种有效、方便的离散方法，但由于电磁场领域的特殊性，这个问题一直没有得到很好的解决。问题的关键在于一方面对复杂的结构，一般的剖分方法难于适用；另一方面，由于剖分的疏密与最终所形成的系数矩阵的存贮量密切相关，因而人们采用了许多方法来减少存储量，如多重网格法，但这些方法的实现较为困难〔6〕。

网格剖分与加密是有限元方法发展的瓶颈之一，采用自适应网格剖分和加密技术相对来说可以较好地解决这一问题。自适应网格剖分根据对场量分布求解后的结果对网格进行增加剖分密度的调整，在网格密集区采用高阶插值函数，以进一步提高精度，在场域分布变化剧烈区域，进行多次加密。

这些年有限元方法的发展日益加快，与其他理论相结合方面也有了新的进展，并取得了相当应用范围的成果，如自适应网格剖分、三维场建模求解、耦合问题、开域问题、高磁性材料及具有磁滞饱和非线性特性介质的处理等，还包括一些尚处于探索阶段的工作，如拟问题、人工智能和专家系统在电磁装置优化设计中的应用、边基有限元法等，这些都使得有限元方法的发展有了质的飞跃。

矩量法将连续方程离散化为代数方程组，既适用于求解微分方程，又适用于求解积分方程。他的求解过程简单，求解步骤统一，应用起来比较方便。然而 77他需要一定的数学技巧，如离散化的程度、基函数与权函数的选取，矩阵求解过程等。另外必须指出的是，矩量法可以达到所需要的精确度，解析部分简单，可计算量很大，即使用高速大容量计算机，计算任务也很繁重。矩量法在天线分析和电磁场散射问题中有比较广泛地应用，已成功用于天线和天线阵的辐射、散射问题、微带和有耗结构分析、非均匀地球上的传播及人体中电磁吸收等。

FDTD用有限差分式替代时域麦克斯韦旋度方程中的微分式，得到关于场分量的有限差分式，针对不同的研究对象，可在不同的坐标系中建模，因而具有这几个优点，容易对复杂媒体建模，通过一次时域分析计算，借助傅里叶变换可以得到整个同带范围内的频率响应；能够实时在现场的空间分布，精确模拟各种辐射体和散射体的辐射特性和散射特性；计算时间短。但是FDTD分析方法由于受到计算机存储容量的限制，其网格空间不能无限制的增加，造成FDTD方法不能适用于较大尺寸，也不能适用于细薄结构的媒质。因为这种细薄结构的最小尺寸比FDTD网格尺寸小很多，若用网格拟和这类细薄结构只能减小网格尺寸，而这必然导致计算机存储容量的加大。因此需要将FDTD与其他技术相结合，目前这种技术正蓬勃发展，如时域积分方程／FDTD方法，FDTD／MOM等。FDTD的应用范围也很广阔，诸如手持机辐射、天线、不同建筑物结构室内的电磁干扰特性研究、微带线等〔7〕。

复射线技术具有物理模型简单、数学处理方便、计算效率高等特点，在复杂目标散射特性分析等应用领域中有重要的研究价值。典型的处理方式是首先将入射平面波离散化为一组波束指向平行的复源点场，通过特定目标情形下的射线追踪、场强计算和叠加各射线场的贡献，可以得到特定观察位置处散射场的高频渐进解。目前已运用复射线分析方法对飞行器天线和天线罩（雷达舱）、（加吸波涂层）翼身结合部和进气道以及涂层的金属平板、角形反射器等典型目标散射特性进行了成功的分析。尽管复射线技术的计算误差可以通过参数调整得到控制，但其本身是一种高频近似计算方法，由于入射波场的离散和只引入鞍点贡献，带来了不可避免的计算误差。总的来说复射线方法在目标电磁散射领域还是具有独特的优势，尤其是对复

杂目标的处理。

5 结 语

电磁学的数值计算方法远远不止以上所举，还有边界元素法、格林函数法等，在具体问题中，应该采用不同的方法，而不应拘泥于这些方法，还可以把这些方法加以综合应用，以达到最佳效果。

电磁学的数值计算是一门计算的艺术，他横跨了多个学科，是数学理论、电磁理论和计算机的有机结合。原则上讲，从直流到光的宽频带范围都属于他的研究范围。为了跟上世界科技发展的需要，应大力进行电磁场的并行计算方法的研究，不断拓广他的应用领域，如生物电磁学、复杂媒质中的电磁正问题和逆问题、医学应用、微波遥感应用、非线性电磁学中的混沌与分叉、微电子学和纳米电子学等。

参考文献

〔1〕 文舸一．计算电磁学的进展与展望〔J〕．电子学报，1995，23（10）：62-69

〔2〕 刘圣民．电磁场的数值方法〔M〕．武汉：华中理工大学出版社，1991．

〔3〕 张成，郑宏兴．小波矩量法求解电磁场积分方程〔J〕．宁夏大学学报（自然科学版），2000，21（1）：76-79

〔4〕 王长清．时域有限差分(FD-TD)法〔J〕．微波学报,1989,(4):8-18

〔5〕 阮颖诤．复射线理论及其应用〔M〕．成都：电子工业出版社，1991．

〔6〕 方静，汪文秉．有限元法和矩量法结合分析背腔天线的辐射特性〔J〕．微波学报，2000，16(2):139-143

〔7〕 杨永侠，王翠玲．电磁场的FDTD分析方法〔J〕．现代电子技术,2001,(11):73-74

〔8〕 洪伟．计算电磁学研究进展〔J〕．东南大学学RB （自然科学版），2002，32（3）：335-339

〔9〕 王长清，祝西里．电磁场计算中的时域有限差分法〔M〕．北京：北京大学出版社，1994．

〔10〕 楼仁海，符果行，袁敬闳．电磁理论〔M〕．成都：电子科技大学出版社，1996．

地震面波频散是研究地球内部结构的有力工具。当前这类反演方法常采用两步反演来实现，首先作区域化反演得到相（群）速度分布，然后作S波速度反演得到速度分布的三维图像。

在本项研究中，为了获得研究区高精度面波层析成像结果，在建立大尺度模型时是按照2°×2°的网度进行，而对于小尺度的精细模型则选用1 °×1 °的网度；所采集信号的周期范围也从8 s扩大到400 s；反演成像时，应用改进后的高精度面波频散反演方法。为了保证反演结果稳定可靠，需要较充足的高质量地震面波资料。所以，我们经过对数字地震记录进行处理筛选得到周期在8～400 s之间2500多条高质量的面波频散曲线，这就为我们进行三维结构面波频散反演奠定了坚实的基础。

面波速度反演一般有两种方法：网格法和非网格法。网格法是把研究区域网格化，每一网格中的速度是均匀的，由实际测得的频散数据来反演每一网格的相速度或群速度。路径的积分由网格分段进行，进行射线路径追踪。宋仲和等（1992）采用此方法反演中国大陆及其相邻海域的面波速度结构及地壳上地幔三维S波速度结构。非网格法是把研究区域的速度函数用一组函数展开。非网格法又分为两种：一种是选择一组先验函数，一般选择一组正交函数；另一种是根据路径的分布，由具体路径来确定一组函数。主要有Taran-tola方法和Yanovskaya-Ditmar方法。Tarantola方法是非线性反演方法，它需要选择模型的先验协方差。

张禹慎和Toshiro Tanimoto（1991），张禹慎和马石庄（1997）利用球谐分析方法研究了全球范围（5°×5°），得到周期在85～250 s的Rayleigh波和Love波的相速度变化。在区域性地震面波层析成像方面，Feng Chi-Chi和邓大量（1983），庄真等（1987）则利用适配滤波频时分析和网格反演方法，分别研究了欧亚大陆和太平洋地区的三维速度结构。在目前还没有一个完整的理论能够计算横向非均匀模型的理论频散的情况下，网格反演方法比较成功地解决了从混合路径频散提取较小单元纯路径频散问题。

5221　面波频散层析成像原理及方法

由于面波具有频散的性质，每个单色波都以它自己的速度传播，而且它主要沿地球表面进行传播。对于瑞利波基阶振型而言，它的每一个频率的波对其波长约1/3深度内的结构很敏感（Knopoff,1972）。因而，利用面波的这些特性，只要在震源和台站间的速度结构没有很强的横向变化，那么，面波就可以被用来研究这两点间的深处的平均结构信息（特别是海洋和体波难以达到的地区）。由此，随着宽频带数字地震仪的发展，面波层析成像也就逐步发展起来。宽频带数据使得速度结构反演分辨率越来越高。

在面波层析成像的传统方法中，一般假定面波绕地球沿大圆路径传播。传统的大圆近似是Evernden（1953,1954）、Mc Garr（1969）、Capon（1970,1971）用LASA和NOR-SAR台阵观测20～30 s周期范围内异常海洋面波的结果后提出的。而根据费马原理，在几何光学近似和一阶扰动理论情况下，该假定是正确的（Nolet,1989）。依据此原理进行层析成像时，区域化是把所研究的每个区域当作是横向均匀的，每个区域具有区域纯路径群速度。沿着一条路径观测到的总相移，作为沿该大圆路径在各区域相移的总和，而这些区域相移已由相应纯路径相对长度加权了。这是一个用区域相速度的“相位积分近似”的大圆路径有限长度的表达式。

在弱横向不均匀的情况下，大量可测的扰动相对于横向均匀介质的情况可用线性关系式的适当精度来表示，该关系式是由速度扰动幂级数中的一阶项给出。则沿着射线L的两点A1和A2之间的相位变化为

中国华北地区岩石圈三维结构及演化

式中，C（s）是相速度。

利用费马原理，忽略掉由于远离大圆的路径畸变的变化，相位扰动可能与沿大圆积分的速度扰动有关：

中国华北地区岩石圈三维结构及演化

式中，C0是平均的大圆相速度。该式就是经典面波层析成像的基础（Nolet,1989）。

基于上述理论，针对于我们所要解决的问题，采用两步法研究地壳上地幔三维速度结构。即首先进行网格化反演，求得面波纯频散，然后由所得的纯频散数据反演华北地区地壳上地幔三维速度结构。具体步骤如下：

1）对原始记录进行预处理，测定面波频散曲线。它包括截取瑞利面波，进行频率时域分析等，以得到实际频散数据。采用频率时域分析法测定群速度频散，对实测地震记录经仪器频响特性校正得到时间信号X（t），在对X（t）进行Fourier变换得到谱函数X（f），以不同的频率点为中心乘以窗函数，作Fourier逆变换求出该频率或周期Tk（k=1，…，M）、ti（i=1，…，N）的瞬时谱振幅Aik。以到时ti（或换算为群速度Ui）为纵坐标，周期Tk为横坐标，绘制矩阵ANM的等值线图，Aik的最大值连线即为所求的群速度频散曲线。

2）对研究区域进行网格划分，计算理论群速度。即采用网格化纯路径频散反演方法得到各块的纯频散数据。

3）为了有效得到S波速度分布，必须选择一个合理的初始模型，本研究利用已有的大地测量、地质、地球物理资料，改进地球圈层结构模型，对每一块建立地壳上地幔的三维速度结构模型作为反演计算的初始模型，使得我们进行该区面波频散反演计算时的S波速度值得到很好的约束。

4）由得到的每块纯频散数据作为初始数据，并利用LSQR及单纯形替换法、SVD方法等广义线形反演方法，对每一块的纯频散曲线进行反演求取该单元的S波速度、P波速度以及密度等随深度的分布。

5）对所得中国及邻区地壳上地幔的三维速度结构进行横向、纵向成像。分析研究该区域地壳上地幔三维速度结构、构造演化及其动力学性质。

为了研究地壳上地幔三维速度结构，我们首先采用网格化纯路径频散反演得到各网格单元的纯频散数据。纯路径频散适于研究大范围内某种分格尺度上的分区平均构造及横向不均匀性的有效方法。基于纯路径频散，首先对研究区域剖分为n个网格单元；并把每个网格看作是横向均匀的，且具有区域纯路径群速度值。然后求取大园路径或进行射线追踪；对某一周期T，根据第j个网格所假定的群速度Uj（T）计算出第i条路径（i=1,2，…，m; m为反演区域内的总射线条数）的理论群速度υi（T）。根据费马原理，在几何光学近似和一阶扰动理论情况下，可以假定面波绕地球沿大圆路径传播。因而，我们首先根据震源和台站位置，求出面波传播的大圆路径，并记下它在所通过的网格中的路径长度。如对于第i条路径，我们可求出它在第j个网格单元内的长度dij以及大圆路径总长度Di。根据上述理论，假定该大圆路径的总相移为各均匀网格单元相移之和，而网格边界上无相移。所以第i条路径的理论群速度υi（T）为

中国华北地区岩石圈三维结构及演化

如果第i条路径实测的群速度为ui（T），而待求的第j（j=1,2，…，n）个网格单元内的区域纯路径群速度值为uj（T），同理可得

中国华北地区岩石圈三维结构及演化

用（54）式减去（53）式，可得出第i条路径实测群速度与理论群速度之差（即群速度的扰动）为

中国华北地区岩石圈三维结构及演化

写成矩阵形式如下：

中国华北地区岩石圈三维结构及演化

其中，b是走时残差向量，A是网格单元内射线长度dij与大圆路径总长度Di比值的大型稀疏矩阵，X是待求的理论群速度与实测群速度倒数之差向量。求解矩阵方程（56）式，对于面波层析成像来说，矩阵A是一个大型稀疏矩阵。这一矩阵通常是奇异的，其逆算子往往不存在，因而必须采用迭代方法，即求其广义逆。通常广义逆有两种方法。一类是最小二乘解意义下的广义逆，即求误差向量L2的最小范数。奇异值分解法（SVD），共轭梯度法（CG），及最小二乘共轭梯度法（LSQR）都属这类方法；另一类寻求Moore—Penrose意义下的广义逆，它不求误差向量L2的极小，而是求模型参量的最小范数解，属于这类方法的有正交化算法ORTH及阻尼最小范数解DMNLS法等。而我们的目的是寻找一种精度高、受数据误差干扰小、收敛快的算法。它应具备：①能充分利用先验信息，对反演参数进行约束，使反演结果更为合理，接近真实模型；②采用数学方法把病态问题转化为良态问题求解；③能充分利用大型稀疏矩阵的特点，减少计算机内存和提高计算精度。通过大量实践证明（曹俊兴，1996），SVD法计算精度高，但计算速度慢，且不能利用大型稀疏矩阵的特点；正交化算法ORTH计算精度与SVD相近，但计算速度也较慢；而最小法二乘共轭梯度法（LSQR，见下述）具有收敛快、稳定性好的特点，以及易于用阻尼因子控制其反演结果质量，且能充分利用大型稀疏矩阵的特点，因此，它被认为是最适合于求解矩阵方程（56）式的算法。

图58 网格化纯路径频散反演方法流程图

根据上述原理，我们设计了一套比较成熟的计算软件（图58），并在实际中得到应用。在编制程序过程中，根据大型稀疏矩阵的特点，对（56）式中的非零元素采取了按行存储格式，这样不但占用内存少，而且也减少了计算量。地震射线的分布不均匀以及台站的影响，使得反演结果出现偏差，因而在程序实现中，我们对射线进行了归并和剔除，对台站附近的单元数据用其周围的单元数据内插得到。这样虽然降低了反演结果的分辨率，但避免了假异常的出现。

LSQR即最小二乘QR 分解（Least Square or Decomposition），它将Lanczos迭代方法结合到共轭梯度（CG）方法中并类似于CG的一种算法，它直接求解Ax=b，而不同于CG的求解 。它也属于迭代算法，在迭代求解过程中使残差范数 单调减少的同时也使 单调减少，因而得到期望解。在求解过程中只涉及非零元素，减少了存储空间，提高运算速度，所以该方法特别适于求解系数为大型稀疏矩阵的方程组，在处理实际资料时可加阻尼来减小噪音。

LSQR求解大型稀疏方程组的算法如下。Paige和Saunders（1982）将Lanczos迭代方法（Lanczos,1950）结合到共轭梯度方法中，提出了一个与CG算法相类似的算法，称之为LSQR方法。LSQ R算法与CG算法的主要差别在于前者求解的是Ax=b而后者求解的是 。求解Ax=b时观测数据误差的放大因子是奇异值的倒数，求解 时观测数据误差的放大因子是奇异值平方的倒数，因此，对于病态问题，LSQ R 算法较CG算法的效果要好，但是，当数据误差较大时LSQR 算法仍会发散。为提高LSQR 算法的抗噪能力，可对迭代求解过程施加一定的阻尼，构成阻尼LSQR 算法。

5222　地震面波频散及波形反演的原始数据处理

为了利用地震面波的频散特性来研究岩石圈三维结构，搜集了90°～140°E,10°～42°N范围内约600个地震事件，并从中挑选出200个事件。这些地震事件起自1982年，终止于2000年，震级大部分都在50～70之间，震源深度小于100km。利用了GDSN、地震台网中部分地震台站的数字化地震记录。这些台站都有三分量的长周期数字地震仪，其采样速率均为1个样点/s。

地震波场数字模拟在地震资料采集布置、处理与解释中具有重要的地位，是地震勘探中的一个有力工具。有限单元法（华东石油学院，1987）和有限差分法（Kelly，et al，1985）等是复杂构造条件下进行地震波场数字模拟的常用方法。它们能够给出地震波动传播的比较详细的过程，记录的波场信息丰富，包括首波、绕射波、折射波、面波、转换波等，但它们占用计算机内存大，耗费计算机时，对模型的大小有一定的限制。常用的射线追踪方法虽具有速度快、节省内存等特点，但它主要反映地震波运动学特点，而不能较好地表现地震波动力学特点，而且存在阴影区。

为克服普通射线法和波动方程方法的某些局限，Cerveny等人相继发展了一种将波动方程和射线理论相结合的方法，被称做高斯射线束法（Cerveny et al，1982，1983，1984，1985；Klimes，1984；Weber，1988；Hill，1990）。该方法同时考虑了波的运动学和动力学特征，适用于复杂的非均匀介质模型，还能考虑介质的吸收作用（Weber，1988）。该法不需要两点射线追踪的计算，具有速度快、精度较高（可与有限差分的计算结果对比（George et al，1987））的特点，对焦散区、临界区及暗区等奇异区域都具有较好的效果。

在国外，有关高斯射线束方法的文章较多，不少研究者对于高斯射线束法合成记录的某些参数选取（Muler，1984；Weber，1988）、合成记录的精度（Weber，1988；George et al，1987）进行了比较深入的研究，对于高斯射线束合成记录的运动学和动力学射线追踪也有不少文章进行了介绍（Cerveny et al，1982；Muler，1984）。Hill（1990）还提出了高斯射线束偏移的方法。尽管这些方法还有一些问题（如参数选取、射线密度等）值得研究，还需要不断完善和更新，但总的看来，它已经形成了地震模型正演的一个方向。

本节将对二维高斯射线束的原理，射线的运动学和动力学追踪，合成记录的实现作较为详细的介绍，文献［6］介绍了计算任意均匀各向同性介质P波、P-SV波地震记录。并导出了利用爆炸反射面原理计算合成叠加剖面的高斯射线束法动力学射线追踪边界条件，及相应的软件，并介绍了大量的研究成果。

431 二维高斯射线束的表达式

在这里，只给出P波位移的高斯射线束表达式，P-SV波有同样的形式。

所谓高斯射线束就是指d性动力学方程集中于射线附近的高频渐近时间调和解，据Cerveny等人（Cerveny et al，1982）的文章，该解在频率域具有如下表达形式

地球物理数据处理教程

式（431）为P波位移的表达式；UP（s，n，ω，t）便是P波位移；ω为圆频率；t为时间参量。

在（431）式中，τ（s）= （s）ds表示波沿中心射线的传播时间，（s，n）为沿某一中心射线Ω的X-Z平面上某点的射线中心坐标，s为从震源沿中心射线到该点的距离，n为垂直于中心射线的距离，如图44所示：

α（s）为波沿中心射线的传播速度；

K（s）=α（s）Re［p（s）/q（s）］为射线的相前曲率；

L（s）= 为射线的有效半宽度；

p（s）、q（s）为随着波沿射线的传播而变化的两复值函数，它们满足：

q，s=αp p，s=-α-2α，nnq （432）

式中q，s= ，同样α，nn= 。

在（431）式中还有一个函数A（s），它是沿中心射线的P波位移幅值。对于层状介质，假定一条射线从炮点S0出发，经过N个界面后反射，回到地面接收点R，则该振幅的表达式为：

地球物理数据处理教程

式中：Ri为界面的反射或透射系数；Q为某界面的反射或透射点；σ为介质的密度；符号“-”表示在射线通过界面时生成射线（反射或透射射线）一侧的量值；αi、βi分别是入射射线和生成射线与局部坐标X轴正方向（界面切线方向）的夹角，如图45所示。

图44 射线中心坐标系示意图

图45 波在界面上的反射与透射示意图

从（431）式中可以看出，位移UP的幅值随着距中心射线的距离n的增大而呈指数衰减，从垂直于中心射线的断面来看，其振幅的分布是钟形的，即呈高斯型分布，这就是高斯射线束名称的由来，L（s）便是一个依赖于频率的射线束之半宽度，如图46所示。从（431）式还可以看出，P波位移的相位除与s有关外，还与n有关，在中心射线上，延时为τ（s）= （s）ds，离开中心射线，延时便与n和K（s）有关，从物理的观点来看，K（s）表示了射线束波前面的曲率。

图46 高斯射线束示意图

从（431）式我们还可以得出沿射线的p、q还应当满足下列条件，使得沿整条射线的高斯射线束处处正则（具有有限振幅）以及保证高斯射线束的集中性（射线有效半宽度是有限的实数）。

（1）q（s）≠0，使得A（s）→/∞ （434）

地球物理数据处理教程

有关二维高斯射线束的表达式（431）的较为详细的推导过程，请见附录A。

432 运动学射线追踪与动力学射线追踪

高斯射线束是波动方程沿射线附近的高频渐近解，它依赖于中心射线，即依赖波的传播路径，为此需要作运动学射线追踪，也即一般简称的射线追踪，而动力学射线追踪则是附加在运动学射线追踪之上的求解高斯射线束表达式中p、q函数的过程，即求解沿射线波的某些动力学特征的传播过程，正是由于动力学射线追踪，才使得高斯射线束具有不同于一般射线法正演的特点。

运动学射线追踪可以采用任何有效的射线追踪方法，在这里，我们给出直接求解程函方程的射线追踪方法。

程函方程为：

地球物理数据处理教程

式中：τ为旅行时；α为波速。

用特征曲线法可得到相应的常微分方程组：

地球物理数据处理教程

式中：Xi（i=1，2，3）表示直角坐标系下的坐标分量，X1=x，X2=y，X3=z；Pi（i=1，2，3）分别为慢度矢量 （cosα，cosβ，cosγ为方向余弦）在直角坐标系下的各分量（Px=P1，Py=P2，Pz=P3）。

在二维介质中，速度α与y轴无关，并假定射线在初始时刻τ=τ0的初始方向仅限于XZ平面内，即Py|τ=τ0=0，式（437）、（438）就减为四个方程

地球物理数据处理教程

现引入角度i、j，如图47所示，i表示射线切线与Z轴的夹角，j表示射线切线在XY平面上投影与X轴方向的夹角，即射线的方向角，有

图47 角度i、j的定义

Px=α-1sini·cosj

Py=α-1sini·sinj （4310）

Pz=α-1cosi

由式（437）、（438）和（4310）有：

地球物理数据处理教程

对于二维情形，上式成为：

地球物理数据处理教程

式（4311）或（4312）是一阶的微分方程组，可用龙格-库塔法求解。由于数字积分是按时间步进行的，到某一特定时间步时，射线可能跨过某速度界面而不是正好在界面上，这时需要求射线与界面的交点，如果积分步长较小，可以将这一段射线当作直线而求与界面的交点。如果积分步长较大，则应考虑射线弯曲的影响，如图48所示。设在某一时间步，计算出的射线段为 ，B点已跨过界面Σ，为比较精确地求 与界面Σ的交点，作过A点、B点射线的切线（射线在A、B点的方向由i角确定）与界面Σ交于A′、B′。用另一积分步长Δt′=Δt· 重新计算射线段 ，直到B点靠近T，就一般而言，只需两三次反复即可求得交点T。射线通过界面应满足反射或透射条件（Snell定律）。由此便求得了射线的运动轨迹，亦即求出了某入射角入射的一条中心射线。

图48 射线与界面交点的求取

有了中心射线，便可以作射线的动力学追踪，即求p、q。函p（s）、q（s）在高斯射线束中起着非常重要的作用，它们决定了高斯射线线束的分布状态，也表征了高频地震波场沿射线传播的动力学特征。

利用ds=αdτ（τ（s）为射线旅行时，α为射线速度）将式（432）写成对旅行时的微分方程有：

q，τ=α2p p，τ=-α-1α，nnq （4313）

再将其变回到直角坐标系，由

地球物理数据处理教程

得到式（4313）在直角坐标系下的表达式：

q，τ=α2p p，τ=-α-1αzzq （4315）

其中 αzz=α，xxcos2i-2α，xzcosi·sini+α，zzsin2i，i 为射线之切线与Z 轴之夹角。

选取两组互相独立的初始条件

地球物理数据处理教程

加上边界条件（Cerveny and Psencik，1979）：

地球物理数据处理教程

同样由龙格-库塔法我们可解出线性独立的两组解

地球物理数据处理教程

从而可以求得复数解

p（x，z）=εp1（x，z）+p2（x，z）q（x，z）=εq1（x，z）+q2（x，z） （4319）

其中ε为任意复常数，待定。

式（4317）中上标“-”表示界面生成射线一侧的量值，α、β的意义同图45所示，且

地球物理数据处理教程

ε的选取对高斯射线束有重要意义。它决定了射线束的半宽度L（s）和射线束波前面的曲率，即决定了高斯射线束的性质，有不少文章对此进行了探讨，如Muler（1984），Weber（1988）。特别是Weber 1988年的文章对ε的选取作了比较深入的研究，并用理论模型实验说明了ε的选取对合成记录精度的影响，在计算中采用

地球物理数据处理教程

这种选择使得射线束在终点具有最小的半宽度，对计算有利。

由方程（4315）的线性，对用爆炸反射面模型模拟合成自激自收叠加剖面记录时的边界条件，经推导有：

地球物理数据处理教程

其中： cosβ（K2sinα- sinβ）+2DnR+ （ -α，n）c o s2β

该条件同时考虑了下行射线和上行射线穿过边界时的边界效应，使动力学射线追踪与运动学射线追踪一样只考虑单程的传播过程，而其动力学特征又与自激自收的传播特征一致，这样使动力学的射线追踪计算速度提高了约一倍。

433 高斯射线束合成地震记录

到此为止，已经得到沿某中心射线附近的d性波动方程高频近似P波位移分量在频率域的表达式——高斯射线束，把从震源展开的各高斯射线束在接收点处叠加而得到频率域的P波位移。因为d性波动方程是线性的，所以按如下形式叠加的表达式也应当近似地满足波动方程：

地球物理数据处理教程

式中：uφ表示初始入射角为φ的高斯束的位移（见式（431））；（s，n）表示与射线有关的射线中心坐标系中接收点R的坐标；Φ（φ）为与初始入射角有关的权函数。

如果式（4322）中Φ（φ）已知，则可以利用式（4322）求出介质中任意一点的波场。Φ（φ）的确定方法是将式（4322）的渐近值与二维线源波动方程精确解的渐近值比较而得出。根据Cerveny等（Cerveny et al，1982）有：

地球物理数据处理教程

还有研究者提出更为复杂的Φ（φ）或Φ（φ，ω）（与ω 有关）的形式（Muller，1984）不在此列出。合成记录的获得有多种方法，用高斯波包法求取，引入震源函数f（t），设f（t）是可积的，其频谱为F（ω）：

地球物理数据处理教程

f（t）应当为高频函数，以保证高斯射线束法的高频近似性。

用u（R，t）表示时间域中的波场，R为接收点，利用傅立叶变换，得到：

地球物理数据处理教程

uφ（R，ω）是入射角为φ的高斯射线束：

地球物理数据处理教程

式中Q为中心射线到达地面的位置，如图46所示：

A（R）≈A（Q）

地球物理数据处理教程

所以：

uφ（R，ω）=A（Q）eiωθ-ωG·e-iωt　（4327）

地球物理数据处理教程

从而式（4325）成为：

地球物理数据处理教程

假定上式中A与ω无关，交换式中的积分次序，并令

地球物理数据处理教程

则式（4329）可写为：

地球物理数据处理教程

式中已将φ的积分限改为φ0至φN了。φ0至φN应包括对检波点R有贡献的所有射线束。函数g叫做波包，它沿射线传播。由于f（t）是高频函数，即

F（ω）=0 0≤ω≤ω0，ω0是某一高频 （4332）

所以g也是高频的波包。

将式（4331）写成离散形式有：

地球物理数据处理教程

Δφ为中心射线入射角的间隔。

式（4333）表示了如下的物理意义：波包g（R，φ）从震源以入射角φ发出，这些波包沿射线传播且与射线紧密相连，它们随着射线的传播连续改变其特性，介质中任意一点的波场是这些波包在该点的叠加。

波包的计算即可以用傅氏变换法，也可以用褶积法求取，下面用一种具有解析形式计算波包的子波函数，该子波函数由下式给出：

f（t）=exp［-（2πfmt/γ）2］·cos（2πfmt+v） （4334）

式中的参数fm，γ，v可根据需要选择，f（t）对应具有高斯包络的谐波载体，其中γ控制包络的宽度，fm为其主频率，该子波一般称作Gabor子波，也叫Puzyrev子波或高斯包络子波。利用该子波形式，最终可得到计算波包的近似解析表达式（Cerveny，1983）：

地球物理数据处理教程

式中f=fm·（1-4πfmG/γ2）。

上式是一个近似公式，适用于 的情形。从式中还可以看出高斯波包的主频为f，并且在时间和空间都具有高斯包络的特点。利用式（4335）可大大提高运算速度。

综上所述，用高斯射线束计算波场时，可分为三步：

（1）作射线追踪，通常的做法是将程函方程化为介质中射线轨迹的微分方程组，然后用龙格-库塔法求足够密的射线轨迹；

（2）作动力学射线追踪，沿射线求高斯射线束的解，即对多条射线求解微分方程组得到p和q的值。

（3）对检波点附近的高斯射线束的贡献进行加权叠加，求得波场值。

434 模型试验及应用实例

图49为高斯射线束法与有限差分法合成记录对比的一个角点模型。图中给出了炮点位置及射线路径，从射线分布看，在X坐标0～2000 m左右和4000～6000 m两个区域无射线穿过，为一般的射线法的阴影区。但用高斯射线法得到的记录（图410）在这两个区域却仍存在着一定能量分布的绕射，与用有限差分法得到的记录（图411）相比是基本一致的。

图49 角点模型及射线路径

图410 角点模型的高斯射线束法合成共炮点记录

图411 角点模型的有限差分法合成炮点记录

（据George et al，1987）

图412 新疆某地区一个实际地质构造断面图

（图中的数值代表各层的速度，单位m/s）

图413 为高斯射线束法计算的叠加剖面图

图412、413为某地实际的构造断面图及用高斯射线束法计算的叠加剖面图。从叠加剖面中可以看出，它很好地反映出了剖面的构造形态。

RMereu BRoy SWinardhi

（Deptof Earth Sciences,University of Western Ontario,London,Ontario,Canada,N6A 5B7）

摘要　过去几年里，在加拿大东部的加拿大地盾开展了一系列长距离地震折射/宽角反射以及近垂直重合反射法的实验研究工作[2]。地震近垂直反射法是一种很有潜力的方法，它使我们能够进一步了解地壳内的详细结构。上述实验的不足之处仍是信噪比较差的问题。这通常是由于地壳内局部侧向不均匀性引起的散射和射线路径延长效应所造成的。在本项研究中，采用了更有效的图像识别法代替传统的CDP映像法。该方法产生了意想不到的效果，大大提高了对地下各主要反射层的分辨能力。地震折射/宽角反射技术利用了组合爆炸，并采用数百个接收器接收。这一方法的实施使我们得以从一个不同的视野去观察和研究地壳。为了优化地下结构的映像，在对数据进行层析时，又发展了适用于侧向非均匀模型参数化的新方法，即三角单元法与延时最小平方反演技术相结合。为了揭示异常的速度结构，从传统的速度解中减去了区域垂直速度梯度，从而获得了具有明显改善的地下结构映像。

关键词　图像识别　层析　地壳　莫霍面　反射　折射　地震映像

1　引言

过去10年里，在加拿大东部穿越加拿大地盾开展了一系列长距离地震折射/宽角反射以及近垂直重合反射法的实验研究工作。主要的研究区域为中央沉积变质带，GrenvilleFront构造带、Sudbury盆地、Midcontinent断裂系等。在1986年进行的Glimpce实验中，对横跨Great lakes的地震测线采用了空气q震源技术，源距为60m。在最近一次，即1992年的Lithoprobe Abitibi-Grenville的陆地地震折射实验中，采用了44个炮点爆破，平均炮间距为30km。每次爆破由一组415仪器记录，记录点之间的距离为10～15km。剖面长度为180～640Km。上述实验构成了一个密集的射线覆盖区，因而使我们能够应用地震层析成像技术更详细地研究我们所感兴趣的构造区的地震波速变化。大多数折射测线都由近垂直反射测线做补充。这些垂直反射测线在陆地上使用了连续震动源，而在湖中使用了空气q震源技术。上述实验的详细情况及分析方法见Green等[1]、Epili和Mereu[2]、Hamilton和Mereu[3]以及Winardhi和Mereu[4]的文章。

2　分析工具

为了有效地分析处理在上述实验中获得的大量数据，特别是地震折射/反射数据，我们设计并研制了一系列的工具。它们是：①能够在喷墨打印机上绘制大型地震剖面的打印机－绘图仪程序；②一套数字滤波器，包括用于三分量数据的极化滤波器和用于CDP数据的图像识别滤波器；③采用自动及人工方法在屏幕上显示并且识别初至波和续至波的计算机交互程序；④在侧向非均匀模型中发射任何类型的射线（P射线、S射线或转换射线）以及多次波射线的交互型射线追踪程序。

3　利用图像识别技术提高信噪比

在大部分开展地壳近垂直地震反射实验的地区，过去都曾经历过某种构造形变，因而地层不再满足均匀层状的条件，即沿侧向和垂向都是非均匀性的。①地震能量在这类介质内传播经过很长射线路径之后，以稍有不同之时间到达各接收站，因而其走时曲线很不光滑。采用某种叠加技术以求提高信噪比的常规处理方法，在许多情形下因非同相叠加，实际上会破坏信号，以致地下界面成像不能令人非常满意。这种现象在一些炮点道集及CDP道集有明显反射体的情况下也会发生。②本项研究着重于以图像识别方法为基础的信号增强技术，该项技术应用于叠前处理流程，用以改善信噪比。Roy和Mereu[5]，根据信号的侧向连续性、振幅、波形、信号频率等特征来鉴别信号和噪声，利用该技术消除CDP道集的噪声，大大降低了叠加剖面上的噪声干扰水平。在信号叠加的方法上也与过去略有不同，即通过在一个小的移动时窗上计算出具有相同极性的信号能量，这样，就连那些非同相性的信号也都考虑在内了。图1a、图1b给出了一条反射剖面的实例。其中图1a是利用传统处理方法得到的映像，而图1b是利用图像识别法得到的改进了的映像。该图所反映的是位于加拿大地盾东南部中央沉积变质带地震剖面中倾斜剪切带的一部分。

4　层析成像分析

利用层析成像技术对大量的地震折射/宽角反射数据进行了分析，从而得到地壳的速度结构映像。具体步骤如下：①建立一个合理的初始模型；②根据该模型用求正问题的方法构制一组理论走时曲线；③用最优化方法对模型进行细调，直至观测走时曲线与理论走时曲线之间的误差达到最小。初始模型的建立是通过对完整的初至Pg波组进行最小二乘法分析后得到的，于是便产生了一个平均的“标准”地壳模型。在我们的实例中，地表波速为613km/s，随深度的增加波速呈线性增大，至地下40Km处波速为710km/s。研究结果表明，在走时曲线上出现的很多局部分散点都是由近地表岩性沿侧向及垂向的变化所致。我们利用延时法将走时数据转换为近地表的速度结构，然后再利用平滑的走时曲线对较深层的结构做层析成像。

图1　横跨加拿大地盾中央沉积变质带的剖面实例

a—传统处理方法的CDP叠加剖面；b—增强地震信号的图像识别技术

有很多模型参数化的方法，例如Cerveny和Psencik[6],Zelt和Smith[7]。每种方法都有其各自的优缺点。在本研究中，我们选择了最新版的“三角单元法”（Mereu,[8]），该方法是将模型剖分为速度梯度为常数的一系列三角单元。之所以选择三角形做模型的基本单元，是因为它对模型结构的侧向变化及倾斜断层非常有效。沿侧向和垂向的速度变化可以表示为

V（x,z）=ax+bz+c

GradV=ai+bj

其中：V（x,z）是位于模型中的任意一点的地震波速度，GradV是每个三角单元内为常数的速度梯度。常数a、b、c是由每个三角形各节点上的速度值确定的。对每个三角形而言，这三个常数，其速度梯度值是不同的。速度场的梯度值为常数，这样就确保了在每个三角单元内追踪的所有的射线路径均为圆弧，因而对整个模型能够很快地追踪射线。射线追踪是通过求解具有三角形边界线性方程的射线弧方程，并在每个交点上应用斯奈尔定律来实现的。三角单元剖分技术的优点是使我们处理的数据量相对减少。当地下结构比较简单时，三角单元可以剖分得大一些；而当地下结构比较复杂时，则需要将三角单元剖分得小一点。因此，模型剖分的方式并不是唯一的。多次的试验结果表明，一般而言，走时数据对模型剖分的不唯一性并不很敏感。其它试验则表明，如果根据炮间距来确定三角单元的数目，如图2所示，则可以获得最佳的单元数。根据时间延迟、自由表面效应、几何扩展效应以及射线在传播路径中的复传递及复反射系数计算出传输函数，就可以构制出合成地震图。计算机自动搜寻模型中所有主要的走时曲线，并对其做数值编码，确定其终点，进而在模型中以适当的角度进行射线追踪，以便计算出波至时间和振幅。地壳走时数据的反演是个非线性反演问题。对于任何自动反演方法，建立模型时都容易出现不稳定解或不规则路径的情况。为了克服这个问题，我们将原始程序修改成为交互式程序。

图2　采用三角单元法对速度模型参数化的实例

该模型中含有一个断裂带，如图3a所示

新程序的特点是用户通过在计算机键盘上简单地敲几个键就可以修改模型的三角形单元及其速度值。对反演的监控是通过即时显示射线路径以及理论走时和观测走时曲线来实现的。在反演过程中，采用了分层剥离法或扩展单元法[9]（White）对数据进行模拟，反演由表层开始，先构制出浅部地壳结构映像，然后逐渐深入到深部地壳，直至莫霍面。

对莫霍面深度上的模型构制，主要利用了观测的PmP和Pn波走时数据，在那些观测不到PmP波或PmP走时曲线的支很短的地区，莫霍面表现为一个厚的过渡带。为了改善地下结构的映像，将速度剖面减去一个平均的标准地壳速度梯度，即可得到异常速度剖面。此外，我们还发展了一种反映模型不确定性的技术，它以射线覆盖的相对密度为基础，并与模拟退火分析相结合。

图3a显示了一个速度沿侧向变化的模型实例，一条大断层将模型中的高速区与低速区隔开。图2表示了如何利用三角单元法对这样一个模型参数化。在该例中利用了10个炮点的数据。图3b为速度异常映像图，该图明显地改进了地下结构映像，突出地表现了断层的存在。图3c为图3a模型的速度不确定剖面。由图可见，最精确的速度值是在射线覆盖最密集处获得的，而在地壳深部和模型两侧边角的射线稀疏处，估算出的波速值准确性较差。

图3　层析成像图

a—速度模型：一条断层将模型划分为一个高速区和一个低速区；b—速度异常映像图：从（a）减去一个标准地壳，该地壳在地表速度为613km/s，莫霍面附近速度为710km/s;c—速度不确定剖面，该不确定值取决于射线覆盖的密度

最近，我们将层析成像法应用于1992年穿越加拿大地盾的Lithoprobe Abitibi-Grenville 4条地震折射剖面分析，取得了较好效果（Winardhi和Mereu,1997）。研究结果中最令人感兴趣的是获得了Grenville Front构造带（Superior和Grenville之间地质单元边界）的构造映像图。该图揭示出一个异常的速度梯度带向东南方向倾伏，穿透地壳，直达莫霍面。

5　结论

现代地震勘探采用了近垂直反射法和宽角反射/折射法。本文介绍了两个非常有效的技术——图像识别法和速度异常层析技术。这些技术的应用改善了地壳的速度结构映像。

致谢　感谢LSPF的KVasudevan和RMaier，以及西安大略大学的BDunn和JBrunet在计算机硬件方面所给予的技术上的帮助，并提供了地震资料处理软件包。本项研究得到Lithoprobe和NSERC的共同资助，资助号为A1793。

（江涛译，许云校）

参考文献

[1]　AGGreen,BMilkereit,ADavidson,CSpencer,DRHutchinson,DRCannon,WFLee,MWAgena,JCBehrendt,WJHinzeCrustal structure of the Grenville Front and adjacent terranesGeology,1988,16:688～792

[2]　DEpili and RFMereuThe Grenville front tectonic zone:Results from the 1986 Great Lakes onshore seismic wide-angle reflection experimentJoural of Geophysical research,1991,96:16335～16348

[3]　DHamilton and RFMereu2-D tomographic imaging across the North American mid-continent rift systemGeophysical Joural international,1993,112:344～358

[4]　SWinardhi and RFMereuCrustal velocity structure of the Superior and Grenville Provinces of the Southeastern Canadian ShieldCanadian Journal of Earth Sciences,In press,1997

[5]　BRoy and RFMereuSignal enhancement using pattern recognition techniques with application to near vertical crustal seismic reflection experimentsGeophysical research letters,1996,23,1849～1852

[6]　VCerveny and IPsencikGaussian beams in two dimensional laterally varying layered structuresRoyal Astronomical Society,Geophysical Journal,1984,78:65～91

[7]　CAZelt and RBSmithSeismic travel-time inversion for 2-D crustal velocity structureGeophysical Joural international,1992,108,16～34

[8]　RFMereuAn interpretation of the seismic-refraction data recorded along profile SJ-6:Morro Bay-Sierra Nevada,CaliforniaIn:interpretation of the SJ-6 seismic reflection/refraction profile,south-central California,USAIn:Proceedings of the 1985 CCSS Workshop on interpretation of seismic Wave Propagation in laterally Heterogeneous Terranes,USGS Open File Report 87-73,AWWalter and WDMooney（Eds）,1987,20～37

[9]　DWhiteTwo dimensional seismic refraction tomographyGeophysical Joural international,1989,97,223～245．

得到ydx-xdy=d(x^2+y^2)时，其实也可看作：xdy-ydx=-d(x^2+y^2)

得通解为：arctan(y/x)=-a-ln(x^2+y^2)

亦即：-arctan(y/x)=a+ln(x^2+y^2)

将(1,0)代入得：a=0，从而有特解为：-arctan(y/x)=ln(x^2+y^2)

关于你的通解涉及了极限问题：arctan(x/y)=c+ln(x^2+y^2)，在(1,0)无意义。

如果取极限，涉及左右极限是否相等。因此对应c又不确定的取值。

但是可有三角知识推一下：

令β=arctan(x/y)，则tanβ=x/y，则tan(π/2-β)=y/x。从而

π/2-β=arctan(y/x)，即：arctan(x/y)+arctan(y/x)=π/2

因此，由arctan(x/y)=c+ln(x^2+y^2)得：-arctan(y/x)=c-π/2+ln(x^2+y^2)

可立即得到：

a=c-π/2，进而得到c=π/2

或者代入初值条件，定出c=π/2

需要指出，上面推导不甚严格，存在漏洞之处就是反正切函数的值域问题。

上面推导都是在假设x/y是正值的情况下作的，故而arctan(x/y)+arctan(y/x)=π/2成立

如果x/y是负值，显然上式明示arctan(y/x)>π/2与反正切函数的值域矛盾。

应改为arctan(x/y)+arctan(y/x)=-π/2

如此就由arctan(x/y)=c+ln(x^2+y^2)得：-arctan(y/x)=c+π/2+ln(x^2+y^2)

代入初值条件，定出c=-π/2。

所以为了避免这个漏洞产生（主要原因是x/y=0，不适用上面两种情况），

尽量选-arctan(y/x)=a+ln(x^2+y^2)作为通解式合适。

同样的道理，如果问(0,1)处的特解，你的答案是优选的。

在GOCAD项目的框架中，已经提出用三角形来模拟极复杂地质界面。这种对三角形面片的选择是基于任何曲面都可以分解成平面或曲线三角形这一事实而决定的，本节中将展示这种分解可以非常高效的用于解决两点射线追踪问题。射线路径的确定基于费马原理——对于给定的发射点、接收点和反射面，要使每条射线的旅行时最小。最小化过程使用曲面三角剖分的单纯形方法迭代来实现。初始射线可以由试射算法、射线偏移算法或弯曲算法来提供。此外，基于GOCAD软件的几何信息数据可以引入动力学信号，为此，用边界曲面定义三维空间的均匀域，并且发展了一种基于有限状态的自动化新算法，用以确定三维空间中任何给定点的对应区域。

有些文献（GFarin,1988;JLGuiziou,AHaas,1988）提出了几种方法用于解决三维两点射线追踪问题。通常这些方法可以给出满意的结果，但当存在复杂非规则地质体时，如正断层、逆断层、盐丘等，它们的速度极慢并且往往不能适应这些不均匀体。

本节中，介绍一种基于GOCAD几何数据结构的新方法。不同物性的交界面（层位）由数据插值得到的三角剖分曲面代表（简写作“T-surface”）。实际工作中，用两步插值过程来构造T-surface：

（1）第一步插值由DSI方法（JLMallet,1989，本书第三章）实现，其目的在于计算三角形顶点的位置使得T-surface与所有有效数据吻合。

（2）第二步插值由Bezier或Gregory方法（GFarin,1988;JAGregory,1980）实现，使得用平滑曲线三角形近似平面三角形。

与基于Bezier，样条或Nurbs的经典方法相比，Gregory的方法允许考虑：

·当前所有有效的不均匀数据（测井数据，地震数据，斜尺数据）；

·某些不精确的数据类型；

·复杂拓扑结构的层位，例如可以考虑一个与盐丘相交的地层。

GOCAD项目的目的不仅在于提供一个有效的复杂地质界面建模工具，它还可以被用于与这些曲面有关的地球物理应用，如射线追踪、偏移、层析……

521　层位的几何建模和地质意义

下面给出适合于射线追踪的界面（层位）表示法要点。在GOCAD项目中，不同介质的交界面用由无序的三角形面元集合构成的界面图形来表示。面元集合的节点为三角形顶点，节点位置由DSI算法得到的。

假设所有层位都包含在一个代表研究区域D的平行面元体中。层｛H1,H2，…，Hm｝将D分割成一个子区域集合｛D1,D2，…，Dm｝对应于独立的均匀介质，为了定义这些区域，我们将界面定向，也就是每层位有两面（正面和负面）。使用GOCAD提供的图形工具，这种定向可以通过交互的方法实现，这样每个区域可以用一个有向界面的子集来定义。例如，一个区域Dj可以由一系列对应于其边界层位的一些面 来定义。对于给定的三维点P，发展了一种基于有限状态，能够自动返回P所属区域ID号的算法。通过对区域的这种定义，可以很容易的定义介质连续性，而这对被穿过的一系列界面所定义的射线的特征研究来说是非常重要的。

注意到，每个层位都至少分割两种介质。对于我们感兴趣d性波传播来说，一个介质平滑变化的区域可以用一个空间函数集合来描述，刻划其d性性质。下面假设每种介质速度为常数。这样在每一区域内射线为一直线并且根据斯涅尔定理在界面处不连续的改变方向。

522　射线追踪问题

设ρ（E,R,Hr）为连接发射点E到接收点R并在层位Hr上反射的一条射线。假设ρ（E,R,Hr）为由对应于地质模型中ρ（E,R,Hr）与层位Hi的交点的n个接触点Ii组成的多边形线：

地质模型计算机辅助设计原理与应用

记σ（E,R,Hr）为对应射线与模型的接触点Ii的（n个）层位Hi系列，称为ρ（E,R,Hr）的“页码”：

地质模型计算机辅助设计原理与应用

根据ρ（E,R,Hr）的定义可知，Hr至少有一次包括于σ（E,R,Hr）中，并且在复杂的地质条件下，层位Hi可以几次出现在σ（E,R,Hr）中。例如，盐丘、透镜体或逆断层等。

对应于射线路径ρ＝ρ（E,R,Hr）的旅行时T（ρ）由下式定义：

地质模型计算机辅助设计原理与应用

这里Vi为射线在包含线段IiIi+1的地质区域Di中的速度，在被线段IiIi+1穿过区域Di（地层）中速度Vi是一常量，并且只要确定IiIi+1的中点所属的区域Di就可以确定这一速度值。

可以看到，T（ρ）是点｛I0，…，Ii，…，In｝的函数，根据费马原理当且仅当ρ（E,R,Hr）为真射线时，这些点对应于T（ρ）的一个局部极值。我们将应用这一性质来求取逼近一个给定初始近似值ρ0（E,R,Hr）的射线ρ（E,R,Hr）。

ρ（E,R,Hr）确定：设ρk（E,R,Hr）为在第k步时ρ（E,R,Hr）的一个近似值，并且让σk（E,R,Hr）为其对应的“页码”：

地质模型计算机辅助设计原理与应用

如果ρk（E,R,Hr）的所有点除Iik外都是固定的，而Iik可以在相应的层位Hik上移动，那么对应于ρk（E,R,Hr）的旅行时可这样表示：

地质模型计算机辅助设计原理与应用

在第（k+1）步上，如果考虑费马原理，可以移动位于Hik的点Iik到 ，并且使T（Iik｜ρk）是最小的，这样得到一个更好的近似值ρk+1＝ρk+1（E,R,Hr）。由上一个近似值ρk（E,R,Hr）导出的射线ρk+1（E,R,Hr）有如下形式：

地质模型计算机辅助设计原理与应用

上面表达式中 的性质将在下节中精确描述。

动态页码。对比文献（VPeireyra,1988;JLGuiziou,AHaas,1988）中提到的一般方法，这里提出的算法允许“页码”σk（E,R,Hr）从第k步到第（k+1）步时改变。这种“页码”的变化由下面的规则来控制：

规则1。层位 一般来说是相同于Hik的，除非 位于Hik的边界，对于最后一种情况建议在下面两种描述中选取其中之一：

（1）如果 是 与另一曲面H的一个连接点，则让 。换句话说，就是曲面 的改变。有赖于基于几何数据的GOCAD结构，可以很容易实现这一点。

（2）如果 不是 与另一曲面H的连接点，那么射线ρ（E,R,Hr）可能穿出了研究区域，这样模型的宽度不足以确定它。在这种情况下，必须放弃ρ（E,R,Hr）的考虑。

规则2。如果新的射线ρk+1（E,R,Hr）与并没有进入页码σk+1（E,R,Hr）的新的层位相交，则有必要在页码σk+1（E,R,Hr）中增加这些层位，并且在ρk+1（E,R,Hr）中增加相应的射线与模型的接触点。为了确定这些新的接触点和其对应的层位，需要测试Pk+1（E,R,Hr）中所有的线段Ii,k+1Ii+1,k+1与地质模型中所有层位的相交。这一 *** 作是非常耗时的，这也是为什么GOCAD数据库允许使用基于八叉树（octree）技术快速算法的原因（JLMallet,1990;YHuang 1990）。

规则3。可能发生这种情况，ρk+1（E,R,Hr）正切于属于页码σk+1（E,R,Hr）但不是Hr的层位 。此时，存在不同于E和R的两个射线与模型的接触点Iα，k+1和Iβ，k+1，并有：

地质模型计算机辅助设计原理与应用

在这种情况下，建议：

·从ρk+1中去掉Iα，k+1和Iβ，k+1，

·从σk+1中去掉Hα，k+1和Hβ，k+1

用单纯形方法寻找 :

对于与初始页码σ0（E,R,Hr）相联系的给定的一个初始近似射线路径ρ0（E,R,Hr），用一种迭代算法来确定射线路径ρ（E,R,Hr），ρ0（E,R,Hr）的逼近值在算法的每一步k中，Iik在Iik上被移动到对应于T（Iik｜ρk）的最小值的点 上。这个最小值通过使用“单纯形”算法来确定，这种算法基于用来定义T-surface的Hik的初始三角形的平滑曲线插值。

应用实例。在图514中，给出了一个由上述方法获得的射线追踪的例子。可以看到地质情况是比较复杂的，特别是包括一个与给定地层相交的盐丘。为了获得较清晰的图像，在图中显示了较少的射线。

图514　使GOCAD产生的几何数据进行射线追踪的例子（Philippe Nobil等，1990）

可以看到盐丘切割了一个层面，层面位于盐丘内部的部分被移动

基于Bezier或样条插值的经典CAD软件的目标是交互地模拟较好的曲面，而不能生成符合地质应用中遇到的复杂数据的曲面。因此，基于这些方法的软件只能生成抽象的地质曲面，而不是与真实地质界面对应的曲面。与这些经典方法相反，在GOCAD项目开发的几何工具允许模拟极复杂的地质体并且可以同时有效地考虑所有的数据。另外，这样获得的模型可以方便的用于开发地球物理应用程序。本节给出的射线追踪算法并不要求使用超级计算机，它可以在工作站上运行，这要归功于GOCAD的几何数据库的结构。

以上就是关于(二)河道砂体正演分析全部的内容，包括:(二)河道砂体正演分析、微分方程的应用、华北地区天然地震面波层析成像等相关内容解答，如果想了解更多相关内容，可以关注我们，你们的支持是我们更新的动力！