CN107843922A - 一种基于地震初至波和反射波走时联合的层析成像方法 - Google Patents

一种基于地震初至波和反射波走时联合的层析成像方法 Download PDF

Info

Publication number
CN107843922A
CN107843922A CN201711417328.3A CN201711417328A CN107843922A CN 107843922 A CN107843922 A CN 107843922A CN 201711417328 A CN201711417328 A CN 201711417328A CN 107843922 A CN107843922 A CN 107843922A
Authority
CN
China
Prior art keywords
mrow
msub
grid
point
walking
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN201711417328.3A
Other languages
English (en)
Inventor
童思友
王忠成
王凯
梁丹
刘张聚
彭崯
葛林苡
郑洪浩
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Ocean University of China
Original Assignee
Ocean University of China
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Ocean University of China filed Critical Ocean University of China
Priority to CN201711417328.3A priority Critical patent/CN107843922A/zh
Publication of CN107843922A publication Critical patent/CN107843922A/zh
Pending legal-status Critical Current

Links

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/282—Application of seismic models, synthetic seismograms
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/34—Displaying seismic recordings or visualisation of seismic data or attributes
    • G01V1/345—Visualisation of seismic data or attributes, e.g. in 3D cubes
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/70—Other details related to processing
    • G01V2210/74—Visualisation of seismic data

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明公开了一种基于地震初至波和反射波走时联合的层析成像方法，其特征在于，包括以下步骤：获取初至波和反射波的走时信息；建立初始的速度模型；根据给定的初始的速度模型进行初至波和反射波的射线追踪，获得各炮到各个接收点的计算走时和射线路径；由上述的走时信息和计算走时之差，以及射线路径构成反演方程组，解反演方程组得到更新后的速度模型；重复上述步骤直到所得到的更新后的速度模型满足精度要求。本发明通过加入反射波，层析成像得到的速度模型更加准确，射线覆盖范围更广，不同层位速度边界更清晰。

Description

一种基于地震初至波和反射波走时联合的层析成像方法
技术领域
本发明属于地球物理反演技术领域，具体地说，涉及一种基于地震初至波和反射波走时联合的层析成像方法。
背景技术
地震走时层析成像是通过观测的地震波传播时间来重建地下速度结构和界面形态的一种技术方法。一直以来，利用初至波层析成像来重建浅地层速度结构是地震勘探领域常用的速度反演方法。初至层析成像是根据地震初至波(包括直达波和折射波)的走时信息进行近地表速度反演的方法。但初至层析成像只利用了地震波的初至走时信息，没有利用反射走时信息。众所周知，地球物理反演是一个非线性反问题，只利用一种地震波的走时进行反演往往会导致反演精度不够，因此尽可能的利用更多的地震走时信息进行反演是一种探索方向。由于初至波的主要由直达波和折射波组成，射线穿过的区域有限，尤其是深层的部分仅有中间能得到覆盖，两侧几乎没有射线穿过，这就导致反演时该区域不能得到很好地速度重构，同时会产生地质假象。此外，由于初至层析成像不需要对地质体划分层位，最终反演出的速度结构边界也不够清晰。
发明内容
有鉴于此，本发明针对上述的问题，提供了一种基于地震初至波和反射波走时联合的层析成像方法，通过加入反射波，层析成像得到的速度模型更加准确，射线覆盖范围更广，不同层位速度边界更清晰。
为了解决上述技术问题，本发明公开了一种基于地震初至波和反射波走时联合的层析成像方法，包括以下步骤：
步骤1、获取初至波和反射波的走时信息；
步骤2、建立初始的速度模型；
步骤3、根据步骤2给定的初始的速度模型进行初至波和反射波的射线追踪，获得各震源点到各个接收点的计算走时和射线路径；
步骤4、由步骤1中的走时信息和计算走时之差，以及射线路径构成反演方程组，解反演方程组得到更新后的速度模型；
步骤5、重复进行步骤3和步骤4直到所得到的更新后的速度模型满足精度要求。
进一步地，所述步骤1中的走时信息包括初至走时和反射走时。
进一步地，所述步骤2中的初始的速度模型为：宽为Q米，深为P米，根据不同的速度分为S层，按照每个网格大小划分为dx米*dz米，将该初始的速度模型用一个M/dx*N/dz的离散网格描述。
进一步地，所述步骤3中的根据步骤2给定的初始的速度模型进行初至波的射线追踪，旅行时线性插值方法进行，具体为：
求某点P最小走时的问题作如下描述：假设A和B两点的坐标和走时已知，求经过线段AB上的某点C到达点P的最小走时；设A点为坐标原点，AB长度为L，AC长度为r，P点坐标为(x，y)，根据旅行时线性插值方法的基本假设，C点的走时TC由A和B两点的走时TA和TB线性插值得到，即：
则P点的走时为C点的走时加上CP线段所用的时间：
式中：ΔT＝TB-TA，s为空间内的慢度；对TP求关于r的导数，并令其等于0得：
当时，
当时，
对一阶导数继续求导，得到TP关于r的二阶导数，看出二阶导数在这两个解处大于0，这两个解都是函数TP的极小值，而TP又是关于r的初等函数，在区间0≤r≤L内连续，对比这2个解以及AB两端点的值，得出TP在0≤r≤L区间内的最小值；
旅行时线性插值实施步骤分为向前计算和向后计算两步，向前计算是从炮点开始，由公式(4-1)和公式(4-2)求出空间内所有网格节点的最小走时；向后计算是从接收点开始，由公式(3)确定射线路径上与各个网格边界的交点。
进一步地，所述初至波射线追踪向前计算的步骤具体为：
ⅰ)计算震源所在网格上4个节点的走时；
ⅱ)计算震源所在列网格上各节点的最小走时，包括以网格下边界为插值线段计算上方2个节点的最小走时，以及以网格上边界为插值线段计算下方2个节点的最小走时，以此类推，直到计算到第一行或最后一行；
ⅲ)计算炮点所在列右侧一列网格节点走时，以网格左边界为插值线段计算网格右方2个节点的最小走时，对于重复计算的节点，若计算所得值比原来的小就替代；
iv)计算炮点所在列右侧一列网格节点走时，由下往上，以网格下边边界为插值线段，计算网格上方2个节点的最小走时，如果比原来的值小就替代；
ⅴ)计算炮点所在列右侧一列网格节点走时，由上往下，以网格上边界为插值线段计算网格下边界2个节点的最小走时，如果比原来的值小就替代；
ⅵ)重复步骤ⅲ-ⅴ，一直到最右边的一列网格；
同理，计算得到炮点左侧网格节点的最小走时。
进一步地，所述初至波射线追踪向后计算的步骤具体为：
ⅰ)根据公式t＝t向前+sd计算经接收点所在网格上4个节点到达接收点的走时，找出其中的最小走时，其中t向前为向前计算得到的网格上节点的走时，s为网格的慢度，d为接收点到节点距离；
ⅱ)根据步骤ⅰ)得到的网格节点，分别以该节点相邻的两条线段为插值线段，用公式(3)、公式(4-1)和公式公式(4-2)计算到达接收点的最小走时以及对应的插值点；
ⅲ)取步骤ⅱ)中求得的两个最小值中较小的一个对应的插值点，作为新的接收点，重复步骤ⅰ)和步骤ⅱ)，得到下一个最小走时对应的插值点；
iv)重复步骤ⅲ)，直到求得的插值点到达炮点所在网格的边界；
ⅴ)从接收点开始，依次连接最小走时所对应的插值点，最后到达震源点，得到了由震源点到接收点最小走时的路径。
进一步地，所述反射波射线追踪向前计算的步骤具体为：
1)先进行下行波走时场的计算，即计算获取由震源点到达反射界面之上每一个网格节点的最小走时，该过程和权利要求5所述初至波射线追踪向前计算部分完全相同。
2)再进行上行波走时场的计算，即计算射线由震源点经过反射界面反射之后到达每一个网格节点的最小走时，具体步骤为：
ⅰ)反射界面所在的节点上的走时取下行走时场所对应的走时；
ⅱ)对于反射界面之上那一行的网格，以网格下边界为插值线段，计算网格上边界2个节点的走时，对于重复计算的节点，如果走时比原来的小就替代；
ⅲ)对于反射界面之上那一行的网格，以网格左边界为插值线段，计算网格右边界2个节点的走时，如果得到的走时比原来的小就替代；
iv)对于反射界面之上那一行的网格，以网格右边界为插值线段，计算网格左边界2个节点的走时，如果得到的走时比原来的小就替代；
ⅴ)现在考虑反射界面之上第二行的网格，重复进行步骤ⅱ)到步骤iv)，如此不断往上推，一直到第一行的网格；
vi)最终可获得上行波的走时场；
所述反射波射线追踪向后计算的步骤具体为：
ⅰ)根据公式t＝t向前+sd计算经接收点所在网格上4个节点到达接收点的走时，找出其中最小的，其中t向前为上行波走时场网格节点的走时，s为网格的慢度，d为接收点到节点距离；
ⅱ)根据步骤ⅰ)得到的网格节点，分别以该节点相邻的两条线段为插值线段，用公式(3)和公式(4)计算到达接收点的最小走时以及对应的插值点；
ⅲ)取步骤ⅱ)中求得的两个最小值中较小的一个对应的插值点，作为新的接收点，重复步骤ⅰ)和步骤ⅱ)，得到下一个最小走时对应的插值点；
iv)重复步骤ⅲ)，直到求得的插值点到达反射界面；
v)把到达反射界面上的差值点作为接收点，根据公式t＝t向前+sd计算经接收点所在网格上4个节点到达接收点的走时，找出其中最小的，其中t向前为下行波走时场网格节点的走时，s为网格的慢度，d为接收点到节点距离；
vi)根据步骤v)得到的网格节点，分别以该节点相邻的两条线段为插值线段，用公式(3)和公式(4)计算到达接收点的最小走时以及对应的插值点；
vii)取步骤vi)中求得的两个最小值中较小的一个对应的插值点，重复步骤v)和步骤vi)，得到下一个最小走时对应的插值点；
viii)重复步骤vii)，直到求得的插值点到达炮点所在网格的边界，
x)从接收点开始，依次连接最小走时所对应的插值点，最后到达震源点，这样就得到了由震源点到接收点反射波的路径。
进一步地，所述步骤4中的由实际走时和计算走时之差，以及射线路径构成反演方程组，解反演方程组得到更新后的速度模型具体为：
地球物理反演问题一般为非线性问题，对于非线性程度较低的反演问题可用以下线性方程组求解：
Gm＝d (5)
式中：G为M×N的系数矩阵，m为N×1的模型向量，d为M×1的观测数据。具体到层析成像领域，G表示射线路径矩阵，m表示离散的速度网格构成的向量，d为观测到的实际地震走时；对于每一条射线，即从一个震源点到一个接收点，会产生射线路径矩阵的一行，以及一个观测到的地震走时；对于层析成像问题需要用迭代的方式逐渐求得逼近真实解的模型m，此时，式(5)用以下方式表达：
GΔm＝Δd (6)
式中：G为当前速度模型下通过射线追踪得到的射线路径矩阵，Δd为观测走时与通过当前模型计算得到的走时之差，Δm为模型修正量。
进一步地，采用联合迭代重建法对公式(6)进行求解，具体如下：
1)假设初始模型的第j个网格中的慢度为sj，第k个炮点对应的接收点有个Nk个；
其中，j＝1,2,…,M(M为网格总数)；k＝1，2，…，source_num(source_num为炮点总数)；
2)利用前面的射线追踪方法，得到该炮点的每个接收点的理论走时Pn
式中aij为由射线追踪得到的第n条射线在第j个网格内射线长度；
3)求出该炮点的每个接收点实际拾取走时Tn与理论走时Pn之差Δtn
Δtn＝Tn-Pn，n＝1，2...，Nk (8)
4)设第n条射线在第j个网格内的慢度修正值cnj，则
设修正值cnj正比于第j个网格内射线通过的路径αnj与该射线长Rn之比，即：
式中αn是第n条射线的比例常数，Rn是第n条射线的全长
将式(10)、(11)代入式(9)，并整理简化，可得
其中，cj为第j个网格的累记修正量，在算第一个炮点之前初始化为0。
5)重复步骤2)，3)，4)，计算第k+1个炮点，直到计算完source_num个炮点；
6)求每个网格的平均修正值，假设计算完source_num个炮点后，经过第j个网格的射线总条数为Yj，则每个网格的平均修正值为
7)用平均修正值对第j个网格的慢度值sj进行修正，即：
sj＝sj+cj (15)
在修正后，sj的值需要受下列物理条件的约束，即
smin≤sj≤smax (16)
若sj＜smin，则取sj＝smin，若sj＞smax，则取sj＝smax。这里smin、smax为介质慢度的范围。
进一步地，所述的精度要求如下：本次迭代计算得到的速度模型和上一次得到的速度模型作差，然后取二范数，如果二范数小于某一个阈值则停止迭代；二范数为向量的长度，相邻两次迭代速度模型差的二范数很小，说明这两个速度模型相差不大，模型已经收敛到一定程度，继续迭代下去速度模型也不会有太大的改变，因此停止迭代。
与现有技术相比，本发明可以获得包括以下技术效果：
在传统的初至层析成像的基础上加入了反射波，由于未知数的个数不变，方程组中个方程的个数增加，使得反演的多解性有所降低。同时，由于反射波的加入，射线覆盖范围更大，射线密度也相应变大，因此重构出来的速度更准确，边界更清晰。
当然，实施本发明的任一产品并不一定需要同时达到以上所述的所有技术效果。
附图说明
此处所说明的附图用来提供对本发明的进一步理解，构成本发明的一部分，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：
图1是本发明旅行时线性插值(LTI)原理示意图；
图2是本发明LTI向前计算确定各节点初至波最小走时的实现过程图；其中，(a)代表由震源点出发计算得到的震源点所在网格上4个节点的走时；(b)代表由(a)得到的4个节点走时出发，分别向上和向下计算出震源点所在列网格节点上的走时，计算采用的公式为权利要求书中公式(4)，以下(c)～(f)都采用该公式计算；(c)代表在震源所在列右边各个网格内，以网格左边界为插值线段，计算得到网格右边界2个节点的走时，对于重复计算的节点，如果所得走时比原来的小则替代；(d)代表在震源所在列右边各个网格内，以网格下边界为插值线段，计算网格上边界2个节点的走时，如果比原来的小就替代；(e)代表在震源所在列右边各个网格内，以网格上边界为插值线段，计算网格下边界2个节点的走时，如果比原来的小就替代；(f)代表整个模型空间内所有网格节点上的最小走时；
图3是本发明LTI向后计算确定初至波射线路径的过程图；其中，(a)代表在接收点所在网格的4个节点上寻找一个初至波走时最小的点，并以该点相邻的两条线段为插值线段，分别寻找一个插值点，使得经过该插值点的路径到达接收点的走时最小，然后对比这两条路径所对应接收点的走时，取接收点走时较小的那条路径所对应的插值点作为新的接收点；(b)代表在新的接收点上重复(a)中的过程；(c)代表重复(b)，直到插值点落在震源点所在网格边界；(d)代表依次连接接收点以及上述计算得到的插值点，直到震源点所在的位置，这样就得到了由震源点到接收点走时最小的射线路径；
图4是本发明用来获取各炮点到接收点初至走时和反射走时的真实速度模型图；(a)代表把反射界面所在位置各节点上赋予下行波在对应位置的走时；(b)代表对于反射界面之上那一行的网格，以网格下边界为插值线段，计算网格上边界2个节点的走时，对于重复计算的节点，如果走时比原来的小就替代；(c)代表对于反射界面之上那一行的网格，以网格左边界为插值线段，计算网格右边界2个节点的走时，如果得到的走时比原来的小就替代；(d)代表对于反射界面之上那一行的网格，以网格右边界为插值线段，计算网格左边界2个节点的走时，如果得到的走时比原来的小就替代；(e)代表网格向上移一行，重复进行(b)～(d)的计算，如此不断向上推移，直到网格的第一行；(f)代表获得的整个区域网格节点上的上行波最小走时，即上行波走时场；
图5是本发明LTI向后计算确定初至波射线路径的过程图，(a)代表在接收点所在网格的4个节点上寻找一个上行波走时最小的点，并以该点相邻的两条边界为插值线段，分别寻找一个插值点，使得经过该插值点的路径到达接收点的走时最小，然后对比这两条路径所对应接收点的走时，取接收点走时较小的那条路径所对应的插值点作为新的接收点；(b)代表在新的接收点上重复(a)中的过程；(c)代表重复(b)，直到插值点落在反射界面上；(d)代表以反射界面上的插值点最为新的接收点，在接收点所在网格的4个节点上寻找一个下行波走时最小的点，并以该点相邻的两条边界为插值线段，分别寻找一个插值点，使得经过该插值点的路径到达接收点的走时最小，然后对比这两条路径所对应接收点的走时，取接收点走时较小的那条路径所对应的插值点作为新的接收点；(e)代表重复(d)，直到插值点落在震源点所在网格边界；(f)代表代表依次连接接收点以及上述计算得到的插值点，直到震源点所在的位置，这样就得到了由震源点反射波的射线路径；
图6是本发明用来获取各震源点到接收点初至走时和反射走时的真实速度模型图；
图7是本发明进行反射波射线追踪时由炮点到反射界面的走时图；
图8是本发明进行反射波射线追踪以反射界面为二次震源向上传播的走时图；
图9是本发明第1炮到各接收点的初至波射线路径图；
图10是本发明第11炮第1层的反射波到个接收点的反射波射线路径图；
图11是本发明第11炮第2层的反射波到个接收点的反射波射线路径图；
图12是本发明用来做层析成像的初始速度模型图；
图13是本发明仅用初至波走时层析反演得到的速度模型图；
图14是本发明初至波和反射波走时联合层析成像得到的速度模型图。
具体实施方式
以下将配合实施例来详细说明本发明的实施方式，藉此对本发明如何应用技术手段来解决技术问题并达成技术功效的实现过程能充分理解并据以实施。
本发明公开了一种基于地震初至波和反射波走时联合的层析成像方法，包括以下步骤：
步骤1、获取初至波和反射波的走时信息，走时信息包括初至走时和反射走时，初至走时是指由震源点发出振动信号开始到接收点接收到初至波所用的时间；反射走时是指由震源点发出振动信号开始到接收点接收到反射波所用的时间；具体的获得方法如下：由震源(一般陆上为炸药震源、海上为气枪或电火花震源)激发地震信号，然后由地震仪控制按照一定规则排列的检波器(检波器所在的位置即为接收点)接收并记录振动信号，最终得到.segy格式的地震记录。把这些.segy格式的地震数据读入到地震数据处理软件中(如ProMAX,GeoEast等)，然后由手动拾取即可获得初至走时和反射走时。
步骤2、给定一个初始的速度模型，作为反演的初始模型，该模型要尽可能的准确，同时划分反射层的层位；
所谓的模型要尽可能的准确，是指所建立的模型要尽可能的和真实的地质信息吻合，即要把已知的信息加入到初始模型的建立当中，比如对于海上地震资料，如果有水深数据，可以用这些水深数据作为第一个反射层的位置，同时把模型关于海水的这一部分设置成水速，根据海域的不同一般在1480m/s～1520m/s之间变化；
具体地，初始的速度模型为：宽为M米，深为N米，根据不同的速度将反射层分为S层，按照每个网格大小划分为dx米*dz米，将该初始的速度模型用一个M/dx*N/dz的离散网格描述。
步骤3、根据步骤2给定的初始的速度模型进行初至波和反射波的射线追踪，获得各炮到各个接收点的计算走时和射线路径；
初至波和反射波射线追踪都分为向前计算和向后计算两个步骤。向前计算可获得相应的走时场，对初至波来讲，就是初至波走时场，对于反射波来讲，分为射线到达反射界面之前的下行波走时场，以及射线经反射界面反射之后的上行波走时场。向后计算是从接收点开始，根据走时场逐步追踪出射线路径，最终到达震源点。其中初至波的向后计算是根据初至波走时场计算出一个插值点，使得经过该插值点到达接收点的路径用时最小，然后把该插值点当做接收点，最终插值点会收敛到震源点所在的网格。反射波的向后计算是先根据上行波走时场计算出一个插值点，使得经过该插值点到达接收点的路径用时最小，然后把该插值点当做新的接收点，最终插值点会落在反射界面上，然后以该点作为接受点，根据下行波的走时场计算出一个插值点，使得经过该插值点到达接收点的路径用时最小，最终插值点会落在炮点所在网格的边界上。无论是向前计算获取走时场，还是向后计算获得插值点，采用的方法都是旅行时线性插值方法。
具体地，初至波射线追踪方法，采用的是旅行时线性插值(LTI)方法；下面介绍LTI基本原理。
如图1所示，求某点P最小走时的问题可以作如下描述：假设A和B两点的坐标和走时已知，求经过线段AB上的某点C到达点P的最小走时。其中线段AB称作差值线段，点C称作插值点。设A点为坐标原点，AB长度为L，AC长度为r，P点坐标为(x，y)，根据LTI的基本假设，C点的走时TC可以由A和B两点的走时TA和TB线性插值得到，即：
则P点的走时为C点的走时加上CP线段所用的时间：
式中：ΔT＝TB-TA，s为空间内的慢度。对TP求关于r的导数，并令其等于0可得：
当时，
当时，
对一阶导数继续求导，可以得到TP关于r的二阶导数，可以看出二阶导数在这两个解处大于0，可见这两个解都是函数TP的极小值，而TP又是关于r的初等函数，在区间0≤r≤L内连续，对比这2个解以及AB两端点的值，就可以得出TP在0≤r≤L区间内的最小值。
下面介绍初至波射线追踪的实现步骤。
采用LTI方法进行初至波射线追踪分为向前计算和向后计算两步。向前计算是从炮点开始，由公式(4)求出空间内所有网格节点的最小走时；向后计算是从接收点开始，由公式(3)确定射线路径上与各个网格边界的交点。
初至波射线追踪向前计算的步骤为：
ⅰ)计算震源所在网格上4个节点的走时，如图2(a)所示；
ⅱ)计算震源所在列网格上各节点的最小走时，包括以网格下边界为插值线段计算上方2个节点的最小走时，以及以网格上边界为插值线段计算下方2个节点的最小走时，以此类推，直到计算到第一行或最后一行，如图2(b)所示；
ⅲ)计算震源点所在列右侧一列网格节点走时，以网格左边界为插值线段计算网格右方2个节点的最小走时，对于重复计算的节点，若计算所得值比原来的小就替代，如图2(c)所示；
iv)计算震源点所在列右侧一列网格节点走时，由下往上，以网格下边边界为插值线段，计算网格上方2个节点的最小走时，如果比原来的值小就替代，如图2(d)所示；
ⅴ)计算炮点所在列右侧一列网格节点走时，由上往下，以网格上边界为插值线段计算网格下边界2个节点的最小走时，如果比原来的值小就替代，如图2(e)所示；
ⅵ)重复步骤ⅲ-ⅴ，一直到最右边的一列网格，如图2(f)。
对于炮点左侧网格节点最小走时的计算，可采用和右侧相似的方式得到。
初至波射线追踪向后计算的过程为：
ⅰ)根据公式t＝t向前+sd计算经接收点所在网格上4个节点到达接收点的走时，找出其中最小的，如图3(a)所示，其中t向前为向前计算得到的网格上节点的走时，s为网格的慢度，d为接收点到节点距离；
ⅱ)根据步骤ⅰ)得到的网格节点，分别以该节点相邻的两条线段为插值线段，用公式(3)和公式(4)计算到达接收点的最小走时以及对应的插值点；
ⅲ)取步骤ⅱ)中求得的两个最小值中较小的一个对应的插值点，将步骤ⅱ)中求得的最小值对应的插值点作为新的接收点，重复步骤ⅰ)和步骤ⅱ)，得到下一个最小走时对应的插值点，如图3(b)所示；
iv)重复步骤ⅲ)，直到求得的插值点到达炮点所在网格的边界，如图3(c)所示；
ⅴ)从接收点开始，依次连接最小走时所对应的插值点，最后到达震源点，这样就得到了由震源点到接收点最小走时的路径。
以上介绍了初至波射线追踪的方法，对于反射波射线追踪，可在初至波射线追踪的基础上进行；反射波的射线追踪仍分为向前计算和先后计算两部分。
反射波射线追踪向前计算分为两个步骤，首先计算下行波的走时场(由震源点到反射界面的走时场)，然后上行波的走时场(由反射界面开始到达各个接收点的走时场)。
下行波走时场的计算和初至波射线追踪向前计算部分完全相同，此处不再赘述。下面介绍上行波走时场的计算方法：
ⅰ)把反射界面所在位置各节点上赋予下行波在对应位置的走时，如图4(a)所示；
ⅱ)对于反射界面之上那一行的网格，以网格下边界为插值线段，计算网格上边界2个节点的走时，对于重复计算的节点，如果走时比原来的小就替代，如图4(b)所示；
ⅲ)对于反射界面之上那一行的网格，以网格左边界为插值线段，计算网格右边界2个节点的走时，如果得到的走时比原来的小就替代，如图4(c)所示；
iv)对于反射界面之上那一行的网格，以网格右边界为插值线段，计算网格左边界2个节点的走时，如果得到的走时比原来的小就替代，如图4(d)所示；
ⅴ)把网格上移一行，重复进行步骤ⅱ)到步骤iv)，如图4(e)所示，如此不断往上推，一直到第一行的网格；
vi)最终可获得上行波的走时场，如图4(f)所示。
反射波射线追踪向后计算的过程为：
ⅰ)根据公式t＝t向前+sd计算经接收点所在网格上4个节点到达接收点的走时，找出其中最小的，如图5(a)所示，其中t向前为上行波走时场网格节点的走时，s为网格的慢度，d为接收点到节点距离；
ⅱ)根据步骤ⅰ)得到的网格节点，分别以该节点相邻的两条线段为插值线段，用公式(3)和公式(4)计算到达接收点的最小走时以及对应的插值点；
ⅲ)取步骤ⅱ)中求得的两个最小值中较小的一个对应的插值点，作为新的接收点，重复步骤ⅰ)和步骤ⅱ)，得到下一个最小走时对应的插值点，如图5(b)所示；
iv)重复步骤ⅲ)，直到求得的插值点到达反射界面，如图5(c)所示；
v)把到达反射界面上的差值点作为接收点，根据公式t＝t向前+sd计算经接收点所在网格上4个节点到达接收点的走时，找出其中最小的，如图5(d)所示，其中t向前为下行波走时场网格节点的走时，s为网格的慢度，d为接收点到节点距离；
vi)根据步骤v)得到的网格节点，分别以该节点相邻的两条线段为插值线段，用公式(3)和公式(4)计算到达接收点的最小走时以及对应的插值点；
vii)取步骤vi)中求得的两个最小值中较小的一个对应的插值点，为新的接收点，重复步骤v)和步骤vi)，得到下一个最小走时对应的插值点，如图5(e)所示；
viii)重复步骤vii)，直到求得的插值点到达炮点所在网格的边界，
x)从接收点开始，依次连接最小走时所对应的插值点，最后到达震源点，这样就得到了由震源点到接收点反射波的路径，如图5(f)所示。
下面给出一个具体的应用例：
图6是给定的一个速度模型，该模型4000m宽，1200m深，共分为3层：0m-400m速度为1500m/s，400m-800m速度为2600m/s，800m-1200m速度为3500m/s。网格大小划分为10m*10m，故该模型可用一个120*400的离散网格描述。首先，利用初至波射线追踪的方法，计算从炮点到达反射界面的最小走时，得到初至走时场(震源坐标为(2000，0))，如图7所示。然后以反射界面作为二次震源，计算从反射界面开始到达反射界面之上各个网格节点的走时，得到二次震源走时场，如图8所示。接下来利用LTI向后计算的方法，以反射点为起点，以上行波走时场为依据向后计算射线路径，直至得到与反射界面的交点，此时得到的是由接收点到反射界面的射线路径。再以该交点为起点，以下行波走时场为依据计算射线路径，直到收敛到震源点，此时得到的是反射界面到震源点的射线路径，这两条射线路径一起构成了反射波的射线路径。
步骤4、由实际走时(步骤1中获取初至波和反射波的走时信息或由本发明向前计算部分走时场中提取的走时)和计算走时(根据初始模型进行向前计算，所得到的接收点所在网格节点上的走时)之差，以及射线路径构成反演方程组，解反演方程组得到更新后的速度模型；具体地，
地球物理反演问题一般为非线性问题，对于非线性程度较低的反演问题可用以下线性方程组求解：
Gm＝d (5)
式中：G为M×N的系数矩阵，m为N×1的模型向量，d为M×1的观测数据。具体到层析成像领域，G表示射线路径矩阵，m表示离散的速度网格构成的向量，d为观测到的实际地震走时。对于每一条射线，即从一个震源点到一个接收点，会产生射线路径矩阵的一行，以及一个观测到的地震走时。对于层析成像问题，没有办法事先知道准确的射线路径矩阵G，因此需要用迭代的方式逐渐求得逼近真实解的模型m。此时，式(5)可用以下方式表达：
GΔm＝Δd (6)
式中：G为当前速度模型下通过射线追踪得到的射线路径矩阵，Δd为观测走时与通过当前模型计算得到的走时之差，Δm为模型修正量。
构造完反演方程组之后，接下来要做的事就是如何解该方程组。一般来说，式(6)是一个混定方程组，没有常规意义下的唯一解。本发明采用的是联合迭代重建法(SIRT)进行的反演，该方法具有较高的计算效率和稳定性。下面介绍SIRT的具体实现方法。
1)假设初始模型的第j个网格中的慢度为sj，第k个炮点对应的接收点有个Nk个；
其中，j＝1,2,…,M(M为网格总数)；k＝1，2，…，source_num(source_num为炮点总数)；
2)利用前面的射线追踪方法，得到该炮点的每个接收点的理论走时Pn
式中aij为由射线追踪得到的第n条射线在第j个网格内射线长度；
3)求出该炮点的每个接收点实际拾取走时Tn与理论走时Pn之差Δtn
Δtn＝Tn-Pn,n＝1,2…,Nk (8)
4)设第n条射线在第j个网格内的慢度修正值cnj，则
设修正值cnj正比于第j个网格内射线通过的路径anj与该射线长Rn之比，即：
式中αn是第n条射线的比例常数，Rn是第n条射线的全长
将式(10)、(11)代入式(9)，并整理简化，可得
其中，cj为第j个网格的累记修正量，在算第一个炮点之前初始化为0。
5)重复步骤2)、3)及4)，计算第k+1个炮点，直到计算完source_num个炮点。
6)求每个网格的平均修正值，假设计算完source_num个炮点后，经过第j个网格的射线总条数为Yj，则每个网格的平均修正值为
7)用平均修正值对第j个网格的慢度值sj进行修正，即：
sj＝sj+cj (15)
在修正后，sj的值需要受下列物理条件的约束，即
smin≤sj≤smax (16)
若sj＜smin，则取sj＝smin，若sj＞smax，则取sj＞smax。这里smin、smax为介质慢度的范围。
步骤5、反复进行步骤3和步骤4和直到所得到的速度模型满足精度要求或达到最高迭代次数。是否满足精度要求可用下列方法判断：本次迭代计算得到的速度模型和上一次得到的速度模型作差，然后取二范数，如果二范数小于某一个阈值(比如10～100，实际计算时可根据反演区域的大小，网格的大小等因素变化)，则停止迭代。二范数可通俗的理解为向量的长度，相邻两次迭代速度模型差的二范数很小，说明这两个速度模型相差不大，模型已经收敛到一定程度，继续迭代下去速度模型也不会有太大的改变，因此停止迭代。设置一个最高的迭代次数(如100次)是为程序安全考虑，为防止模型不收敛而导致程序一直计算下去。
实施例1
在图4模型的基础上，在地表面0m，200m，400m，……，3800m处共设置了20个炮点，对于每一个炮点，都有0m，10m，20m，……，4000m共401个接收点与之对应。分别以这20个点为震源点，以400m和800m这两个深度层为反射界面，对每一个接收点进行初至波走时和来自于两个反射层反射波走时进行计算，得到了每一炮对应每一接收点的初至波走时和反射波走时。图9为第1炮(炮点坐标为(0,0))对应各接收点的初至波射线路径，由图可以看出初至波包括直达波以及来自第一层和第二层的折射波。图10为第11炮(炮点坐标为(2000,0))对应各接收点来自第一层的反射波射线路径，图11为第11炮对应各接收点来自第二层的反射波射线路径。由图8和图9可以看出，由于模型相对比较规则，射线路径也比较规则，满足反射定律。
现在以由图6速度模型和以上观测系统得到的走时为观测走时，即对应于野外实际数据的拾取走时，以图12中的速度模型为初始速度模型进行层析反演。其中图12中的速度模型为一个速度为3500m/s的匀速模型。图13是仅用初至波走时进行反演，迭代20次得到的结果，由图可以看出，通过初至层析成像，能够恢复速度模型的轮廓，尤其是对于浅层的速度比较准确。但同时初至层析也有其缺点：速度边界模糊，尤其是第一个反射层的边界很不清楚；对于第二层左下方和右下方，由于射线路径覆盖不到，没有得到很好地反演，因而速度相对不准，带来了反演的假象。图14是用初至波走时和来自两个反射层的反射波走时联合层析成像得到的速度模型，由图可以看出，该速度模型和真实速度模型(图6)相差不大，并且相对于仅用初至波走时层析反演的结果(图13)该模型不同速度边界更清楚、准确，且由于反射波的加入，射线覆盖范围更大，第二层左下角和右下角部分得到覆盖，反演效果明显好于图13。
上述说明示出并描述了发明的若干优选实施例，但如前所述，应当理解发明并非局限于本文所披露的形式，不应看作是对其他实施例的排除，而可用于各种其他组合、修改和环境，并能够在本文所述发明构想范围内，通过上述教导或相关领域的技术或知识进行改动。而本领域人员所进行的改动和变化不脱离发明的精神和范围，则都应在发明所附权利要求的保护范围内。

Claims (10)

1.一种基于地震初至波和反射波走时联合的层析成像方法，其特征在于，包括以下步骤：
步骤1、获取初至波和反射波的走时信息；
步骤2、建立初始的速度模型；
步骤3、根据步骤2给定的初始的速度模型进行初至波和反射波的射线追踪，获得各震源点到各个接收点的计算走时和射线路径；
步骤4、由步骤1中的走时信息和计算走时之差，以及射线路径构成反演方程组，解反演方程组得到更新后的速度模型；
步骤5、重复进行步骤3和步骤4直到所得到的更新后的速度模型满足精度要求。
2.根据权利要求1所述的基于地震初至波和反射波走时联合的层析成像方法，其特征在于，所述步骤1中的走时信息包括初至走时和反射走时。
3.根据权利要求1所述的基于地震初至波和反射波走时联合的层析成像方法，其特征在于，所述步骤2中的初始的速度模型为：宽为Q米，深为P米，根据不同的速度分为S层，按照每个网格大小划分为dx米*dz米，将该初始的速度模型用一个M/dx*N/dz的离散网格描述。
4.根据权利要求1所述的基于地震初至波和反射波走时联合的层析成像方法，其特征在于，所述步骤3中的根据步骤2给定的初始的速度模型进行初至波的射线追踪，旅行时线性插值方法进行，具体为：
求某点P最小走时的问题作如下描述：假设A和B两点的坐标和走时已知，求经过线段AB上的某点C到达点P的最小走时；设A点为坐标原点，AB长度为L，AC长度为r，P点坐标为(x，y)，根据旅行时线性插值方法的基本假设，C点的走时Tc由A和B两点的走时TA和TB线性插值得到，即：
<mrow> <msub> <mi>T</mi> <mi>C</mi> </msub> <mo>=</mo> <msub> <mi>T</mi> <mi>A</mi> </msub> <mo>+</mo> <mfrac> <mi>r</mi> <mi>L</mi> </mfrac> <mrow> <mo>(</mo> <msub> <mi>T</mi> <mi>B</mi> </msub> <mo>-</mo> <msub> <mi>T</mi> <mi>A</mi> </msub> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow>
则P点的走时为C点的走时加上CP线段所用的时间：
<mrow> <mtable> <mtr> <mtd> <mrow> <msub> <mi>T</mi> <mi>P</mi> </msub> <mo>=</mo> <msub> <mi>T</mi> <mi>C</mi> </msub> <mo>+</mo> <mi>s</mi> <msqrt> <mrow> <msup> <mrow> <mo>(</mo> <mi>x</mi> <mo>-</mo> <mi>r</mi> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msup> <mi>y</mi> <mn>2</mn> </msup> </mrow> </msqrt> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mo>=</mo> <msub> <mi>T</mi> <mi>A</mi> </msub> <mo>+</mo> <mfrac> <mi>r</mi> <mi>L</mi> </mfrac> <mi>&amp;Delta;</mi> <mi>T</mi> <mo>+</mo> <mi>s</mi> <msqrt> <mrow> <msup> <mrow> <mo>(</mo> <mi>x</mi> <mo>-</mo> <mi>r</mi> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msup> <mi>y</mi> <mn>2</mn> </msup> </mrow> </msqrt> </mrow> </mtd> </mtr> </mtable> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </mrow>
式中：ΔT＝TB-TA，s为空间内的慢度；对TP求关于r的导数，并令其等于0得：
<mrow> <mi>r</mi> <mo>=</mo> <mi>x</mi> <mo>&amp;PlusMinus;</mo> <mfrac> <mrow> <mi>y</mi> <mi>&amp;Delta;</mi> <mi>T</mi> </mrow> <msqrt> <mrow> <msup> <mi>L</mi> <mn>2</mn> </msup> <msup> <mi>s</mi> <mn>2</mn> </msup> <mo>-</mo> <msup> <mi>&amp;Delta;T</mi> <mn>2</mn> </msup> </mrow> </msqrt> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </mrow>
当时，
<mrow> <msub> <mi>T</mi> <mi>P</mi> </msub> <mo>=</mo> <msub> <mi>T</mi> <mi>A</mi> </msub> <mo>+</mo> <mfrac> <mi>r</mi> <mi>L</mi> </mfrac> <mi>&amp;Delta;</mi> <mi>T</mi> <mo>+</mo> <mfrac> <mrow> <mi>y</mi> <msup> <mrow> <mo>(</mo> <msup> <mi>L</mi> <mn>2</mn> </msup> <msup> <mi>s</mi> <mn>2</mn> </msup> <mo>+</mo> <msup> <mi>&amp;Delta;T</mi> <mn>2</mn> </msup> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> <mrow> <mi>L</mi> <mrow> <mo>(</mo> <msup> <mi>L</mi> <mn>2</mn> </msup> <msup> <mi>s</mi> <mn>2</mn> </msup> <mo>-</mo> <msup> <mi>&amp;Delta;T</mi> <mn>2</mn> </msup> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>4</mn> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow>
当时，
<mrow> <msub> <mi>T</mi> <mi>P</mi> </msub> <mo>=</mo> <msub> <mi>T</mi> <mi>A</mi> </msub> <mo>+</mo> <mfrac> <mi>r</mi> <mi>L</mi> </mfrac> <mi>&amp;Delta;</mi> <mi>T</mi> <mo>+</mo> <mfrac> <mi>y</mi> <mi>L</mi> </mfrac> <msqrt> <mrow> <msup> <mi>L</mi> <mn>2</mn> </msup> <msup> <mi>s</mi> <mn>2</mn> </msup> <mo>-</mo> <msup> <mi>&amp;Delta;T</mi> <mn>2</mn> </msup> </mrow> </msqrt> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>4</mn> <mo>-</mo> <mn>2</mn> <mo>)</mo> </mrow> </mrow>
对一阶导数继续求导，得到TP关于r的二阶导数，看出二阶导数在这两个解处大于0，这两个解都是函数TP的极小值，而TP又是关于r的初等函数，在区间0≤r≤L内连续，对比这2个解以及AB两端点的值，得出TP在0≤r≤L区间内的最小值；
旅行时线性插值实施步骤分为向前计算和向后计算两步，向前计算是从炮点开始，由公式(4-1)和公式(4-2)求出空间内所有网格节点的最小走时；向后计算是从接收点开始，由公式(3)确定射线路径上与各个网格边界的交点。
5.根据权利要求4所述的基于地震初至波和反射波走时联合的层析成像方法，其特征在于，所述初至波射线追踪向前计算的步骤具体为：
ⅰ)计算震源所在网格上4个节点的走时；
ⅱ)计算震源所在列网格上各节点的最小走时，包括以网格下边界为插值线段计算上方2个节点的最小走时，以及以网格上边界为插值线段计算下方2个节点的最小走时，以此类推，直到计算到第一行或最后一行；
ⅲ)计算炮点所在列右侧一列网格节点走时，以网格左边界为插值线段计算网格右方2个节点的最小走时，对于重复计算的节点，若计算所得值比原来的小就替代；
iv)计算炮点所在列右侧一列网格节点走时，由下往上，以网格下边边界为插值线段，计算网格上方2个节点的最小走时，如果比原来的值小就替代；
ⅴ)计算炮点所在列右侧一列网格节点走时，由上往下，以网格上边界为插值线段计算网格下边界2个节点的最小走时，如果比原来的值小就替代；
ⅵ)重复步骤ⅲ-ⅴ，一直到最右边的一列网格；
同理，计算得到炮点左侧网格节点的最小走时。
6.根据权利要求4所述的基于地震初至波和反射波走时联合的层析成像方法，其特征在于，所述初至波射线追踪向后计算的步骤具体为：
ⅰ)根据公式t＝t向前+sd计算经接收点所在网格上4个节点到达接收点的走时，找出其中的最小走时，其中t向前为向前计算得到的网格上节点的走时，s为网格的慢度，d为接收点到节点距离；
ⅱ)根据步骤ⅰ)得到的网格节点，分别以该节点相邻的两条线段为插值线段，用公式(3)、公式(4-1)和公式公式(4-2)计算到达接收点的最小走时以及对应的插值点；
ⅲ)取步骤ⅱ)中求得的两个最小值中较小的一个对应的插值点，作为新的接收点，重复步骤ⅰ)和步骤ⅱ)，得到下一个最小走时对应的插值点；
iv)重复步骤ⅲ)，直到求得的插值点到达炮点所在网格的边界；
ⅴ)从接收点开始，依次连接最小走时所对应的插值点，最后到达震源点，得到了由震源点到接收点最小走时的路径。
7.根据权利要求1所述的基于地震初至波和反射波走时联合的层析成像方法，其特征在于，所述反射波射线追踪向前计算的步骤具体为：
1)先进行下行波走时场的计算，即计算获取由震源点到达反射界面之上每一个网格节点的最小走时，该过程和权利要求5所述初至波射线追踪向前计算部分完全相同。
2)再进行上行波走时场的计算，即计算射线由震源点经过反射界面反射之后到达每一个网格节点的最小走时，具体步骤为：
ⅰ)反射界面所在的节点上的走时取下行走时场所对应的走时；
ⅱ)对于反射界面之上那一行的网格，以网格下边界为插值线段，计算网格上边界2个节点的走时，对于重复计算的节点，如果走时比原来的小就替代；
ⅲ)对于反射界面之上那一行的网格，以网格左边界为插值线段，计算网格右边界2个节点的走时，如果得到的走时比原来的小就替代；
iv)对于反射界面之上那一行的网格，以网格右边界为插值线段，计算网格左边界2个节点的走时，如果得到的走时比原来的小就替代；
ⅴ)现在考虑反射界面之上第二行的网格，重复进行步骤ⅱ)到步骤iv)，如此不断往上推，一直到第一行的网格；
vi)最终可获得上行波的走时场；
所述反射波射线追踪向后计算的步骤具体为：
ⅰ)根据公式t＝t向前+sd计算经接收点所在网格上4个节点到达接收点的走时，找出其中最小的，其中t向前为上行波走时场网格节点的走时，s为网格的慢度，d为接收点到节点距离；
ⅱ)根据步骤ⅰ)得到的网格节点，分别以该节点相邻的两条线段为插值线段，用公式(3)和公式(4)计算到达接收点的最小走时以及对应的插值点；
ⅲ)取步骤ⅱ)中求得的两个最小值中较小的一个对应的插值点，作为新的接收点，重复步骤ⅰ)和步骤ⅱ)，得到下一个最小走时对应的插值点；
iv)重复步骤ⅲ)，直到求得的插值点到达反射界面；
v)把到达反射界面上的差值点作为接收点，根据公式t＝t向前+sd计算经接收点所在网格上4个节点到达接收点的走时，找出其中最小的，其中t向前为下行波走时场网格节点的走时，s为网格的慢度，d为接收点到节点距离；
vi)根据步骤v)得到的网格节点，分别以该节点相邻的两条线段为插值线段，用公式(3)和公式(4)计算到达接收点的最小走时以及对应的插值点；
vii)取步骤vi)中求得的两个最小值中较小的一个对应的插值点，重复步骤v)和步骤vi)，得到下一个最小走时对应的插值点；
viii)重复步骤vii)，直到求得的插值点到达炮点所在网格的边界，
x)从接收点开始，依次连接最小走时所对应的插值点，最后到达震源点，这样就得到了由震源点到接收点反射波的路径。
8.根据权利要求1所述的基于地震初至波和反射波走时联合的层析成像方法，其特征在于，所述步骤4中的由实际走时和计算走时之差，以及射线路径构成反演方程组，解反演方程组得到更新后的速度模型具体为：
地球物理反演问题一般为非线性问题，对于非线性程度较低的反演问题可用以下线性方程组求解：
Gm＝d (5)
式中：G为M×N的系数矩阵，m为N×1的模型向量，d为M×1的观测数据。具体到层析成像领域，G表示射线路径矩阵，m表示离散的速度网格构成的向量，d为观测到的实际地震走时；对于每一条射线，即从一个震源点到一个接收点，会产生射线路径矩阵的一行，以及一个观测到的地震走时；对于层析成像问题需要用迭代的方式逐渐求得逼近真实解的模型m，此时，式(5)用以下方式表达：
GΔm＝Δd (6)
式中：G为当前速度模型下通过射线追踪得到的射线路径矩阵，Δd为观测走时与通过当前模型计算得到的走时之差，Δm为模型修正量。
9.根据权利要求8所述的基于地震初至波和反射波走时联合的层析成像方法，其特征在于，采用联合迭代重建法对公式(6)进行求解，具体如下：
1)假设初始模型的第j个网格中的慢度为sj，第k个炮点对应的接收点有个Nk个；
其中，j＝1,2,…,M(M为网格总数)；k＝1，2，…，source_num(source_num为炮点总数)；
2)利用前面的射线追踪方法，得到该炮点的每个接收点的理论走时Pn
<mrow> <msub> <mi>P</mi> <mi>n</mi> </msub> <mo>=</mo> <msubsup> <mo>&amp;Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>M</mi> </msubsup> <msub> <mi>a</mi> <mrow> <mi>i</mi> <mi>j</mi> </mrow> </msub> <msub> <mi>s</mi> <mi>j</mi> </msub> <mo>,</mo> <mi>n</mi> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mn>2</mn> <mo>,</mo> <mn>...</mn> <mo>,</mo> <msub> <mi>N</mi> <mi>k</mi> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>7</mn> <mo>)</mo> </mrow> </mrow>
式中αij为由射线追踪得到的第n条射线在第j个网格内射线长度；
3)求出该炮点的每个接收点实际拾取走时Tn与理论走时Pn之差Δtn
Δtn＝Tn-Pn,n＝1,2…,Nk (8)
4)设第n条射线在第j个网格内的慢度修正值cnj，则
<mrow> <msubsup> <mo>&amp;Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>M</mi> </msubsup> <msub> <mi>a</mi> <mrow> <mi>n</mi> <mi>j</mi> </mrow> </msub> <msub> <mi>c</mi> <mrow> <mi>n</mi> <mi>j</mi> </mrow> </msub> <mo>=</mo> <msub> <mi>&amp;Delta;t</mi> <mi>n</mi> </msub> <mo>,</mo> <mi>n</mi> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mn>2...</mn> <mo>,</mo> <msub> <mi>N</mi> <mi>k</mi> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>9</mn> <mo>)</mo> </mrow> </mrow>
设修正值cnj正比于第j个网格内射线通过的路径αnj与该射线长Rn之比，即：
<mrow> <msub> <mi>c</mi> <mrow> <mi>n</mi> <mi>j</mi> </mrow> </msub> <mo>=</mo> <msub> <mi>&amp;alpha;</mi> <mi>n</mi> </msub> <mfrac> <msub> <mi>a</mi> <mrow> <mi>n</mi> <mi>j</mi> </mrow> </msub> <msub> <mi>R</mi> <mi>n</mi> </msub> </mfrac> <mo>,</mo> <mi>n</mi> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mn>2...</mn> <mo>,</mo> <msub> <mi>N</mi> <mi>k</mi> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>10</mn> <mo>)</mo> </mrow> </mrow>
式中αn是第n条射线的比例常数，Rn是第n条射线的全长
<mrow> <msub> <mi>R</mi> <mi>n</mi> </msub> <mo>=</mo> <msubsup> <mo>&amp;Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>M</mi> </msubsup> <msub> <mi>a</mi> <mrow> <mi>n</mi> <mi>j</mi> </mrow> </msub> <mo>,</mo> <mi>n</mi> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mn>2...</mn> <mo>,</mo> <msub> <mi>N</mi> <mi>k</mi> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>11</mn> <mo>)</mo> </mrow> </mrow>
将式(10)、(11)代入式(9)，并整理简化，可得
<mrow> <msub> <mi>c</mi> <mrow> <mi>n</mi> <mi>j</mi> </mrow> </msub> <mo>=</mo> <msub> <mi>&amp;Delta;t</mi> <mi>n</mi> </msub> <mfrac> <msub> <mi>a</mi> <mrow> <mi>n</mi> <mi>j</mi> </mrow> </msub> <mrow> <msubsup> <mo>&amp;Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>M</mi> </msubsup> <msup> <mrow> <mo>(</mo> <msub> <mi>a</mi> <mrow> <mi>n</mi> <mi>j</mi> </mrow> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> </mfrac> <mo>,</mo> <mi>n</mi> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mn>2...</mn> <mo>,</mo> <msub> <mi>N</mi> <mi>k</mi> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>12</mn> <mo>)</mo> </mrow> </mrow>
<mrow> <msub> <mi>c</mi> <mi>j</mi> </msub> <mo>=</mo> <msub> <mi>c</mi> <mi>j</mi> </msub> <mo>+</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>n</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mi>k</mi> </msub> </munderover> <msub> <mi>c</mi> <mrow> <mi>n</mi> <mi>j</mi> </mrow> </msub> <mo>,</mo> <mi>n</mi> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mn>2...</mn> <mo>,</mo> <msub> <mi>N</mi> <mi>k</mi> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>13</mn> <mo>)</mo> </mrow> </mrow>
其中，cj为第j个网格的累记修正量，在算第一个炮点之前初始化为0。
5)重复步骤2)，3)，4)，计算第k+1个炮点，直到计算完source_num个炮点；
6)求每个网格的平均修正值，假设计算完source_num个炮点后，经过第j个网格的射线总条数为Yj，则每个网格的平均修正值为
<mrow> <msub> <mi>c</mi> <mi>j</mi> </msub> <mo>=</mo> <mfrac> <msub> <mi>c</mi> <mi>j</mi> </msub> <msub> <mi>Y</mi> <mi>j</mi> </msub> </mfrac> <mo>,</mo> <mi>n</mi> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mn>2...</mn> <mo>,</mo> <msub> <mi>N</mi> <mi>k</mi> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>14</mn> <mo>)</mo> </mrow> </mrow>
7)用平均修正值对第j个网格的慢度值sj进行修正，即：
sj＝sj+cj (15)
在修正后，sj的值需要受下列物理条件的约束，即
smin≤sj≤smax (16)
若sj＜smin，则取sj＝smin，若sj＞smax，则取sj＝smax。这里smin、smax为介质慢度的范围。
10.根据权利要求1所述的基于地震初至波和反射波走时联合的层析成像方法，其特征在于，所述的精度要求如下：本次迭代计算得到的速度模型和上一次得到的速度模型作差，然后取二范数，如果二范数小于某一个阈值则停止迭代；二范数为向量的长度，相邻两次迭代速度模型差的二范数很小，说明这两个速度模型相差不大，模型已经收敛到一定程度，继续迭代下去速度模型也不会有太大的改变，因此停止迭代。
CN201711417328.3A 2017-12-25 2017-12-25 一种基于地震初至波和反射波走时联合的层析成像方法 Pending CN107843922A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201711417328.3A CN107843922A (zh) 2017-12-25 2017-12-25 一种基于地震初至波和反射波走时联合的层析成像方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201711417328.3A CN107843922A (zh) 2017-12-25 2017-12-25 一种基于地震初至波和反射波走时联合的层析成像方法

Publications (1)

Publication Number Publication Date
CN107843922A true CN107843922A (zh) 2018-03-27

Family

ID=61684113

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201711417328.3A Pending CN107843922A (zh) 2017-12-25 2017-12-25 一种基于地震初至波和反射波走时联合的层析成像方法

Country Status (1)

Country Link
CN (1) CN107843922A (zh)

Cited By (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108983289A (zh) * 2018-05-02 2018-12-11 中国石油天然气集团有限公司 一种确定地震波走时的方法及装置
CN109061725A (zh) * 2018-06-22 2018-12-21 中国电建集团贵阳勘测设计研究院有限公司 弹性波数据采集起跳点的自动识别方法及其所用的设备
CN109100792A (zh) * 2018-10-31 2018-12-28 中国石油化工股份有限公司 基于台站与三维地震联合采集资料的速度反演方法
CN110187382A (zh) * 2019-03-05 2019-08-30 中国石油大学(华东) 一种回折波和反射波波动方程旅行时反演方法
CN111638551A (zh) * 2019-03-01 2020-09-08 中国石油天然气集团有限公司 地震初至波走时层析方法及装置
CN112305595A (zh) * 2019-07-24 2021-02-02 中国石油化工股份有限公司 基于折射波分析地质体结构的方法及存储介质
CN112485825A (zh) * 2019-09-11 2021-03-12 中国石油化工股份有限公司 一种基于初至波走时层析的微测井解释方法
CN112596103A (zh) * 2020-11-24 2021-04-02 中国地质科学院地球物理地球化学勘查研究所 射线追踪方法、装置和电子设备
CN113589374A (zh) * 2020-04-30 2021-11-02 中国石油化工股份有限公司 基于射线密度的有效速度提取方法
CN113777654A (zh) * 2021-08-06 2021-12-10 同济大学 一种基于伴随状态法初至波走时层析的海水速度建模方法
CN113791447A (zh) * 2021-10-12 2021-12-14 同济大学 一种反射结构导引的反射波层析反演方法
CN114879249A (zh) * 2022-04-13 2022-08-09 中国海洋大学 基于四面体单元走时扰动插值的地震波前走时计算方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102967189A (zh) * 2012-11-22 2013-03-13 中北大学 爆炸冲击波超压时空场重建方法
US8406081B2 (en) * 2009-09-25 2013-03-26 Landmark Graphics Corporation Seismic imaging systems and methods employing tomographic migration-velocity analysis using common angle image gathers
CN106353793A (zh) * 2015-07-17 2017-01-25 中国石油化工股份有限公司 一种基于走时增量双线性插值射线追踪的井间地震层析反演方法
CN106680870A (zh) * 2017-03-16 2017-05-17 中国海洋大学 高精度地震波走时射线追踪方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8406081B2 (en) * 2009-09-25 2013-03-26 Landmark Graphics Corporation Seismic imaging systems and methods employing tomographic migration-velocity analysis using common angle image gathers
CN102967189A (zh) * 2012-11-22 2013-03-13 中北大学 爆炸冲击波超压时空场重建方法
CN106353793A (zh) * 2015-07-17 2017-01-25 中国石油化工股份有限公司 一种基于走时增量双线性插值射线追踪的井间地震层析反演方法
CN106680870A (zh) * 2017-03-16 2017-05-17 中国海洋大学 高精度地震波走时射线追踪方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
姜少之: "初至波三维地震CT成像", 《中国优秀硕士学位论文全文数据库 基础科学辑》 *

Cited By (18)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108983289A (zh) * 2018-05-02 2018-12-11 中国石油天然气集团有限公司 一种确定地震波走时的方法及装置
CN108983289B (zh) * 2018-05-02 2020-06-09 中国石油天然气集团有限公司 一种确定地震波走时的方法及装置
CN109061725A (zh) * 2018-06-22 2018-12-21 中国电建集团贵阳勘测设计研究院有限公司 弹性波数据采集起跳点的自动识别方法及其所用的设备
CN109061725B (zh) * 2018-06-22 2020-05-05 中国电建集团贵阳勘测设计研究院有限公司 弹性波数据采集起跳点的自动识别方法及其所用的设备
CN109100792A (zh) * 2018-10-31 2018-12-28 中国石油化工股份有限公司 基于台站与三维地震联合采集资料的速度反演方法
WO2020087767A1 (zh) * 2018-10-31 2020-05-07 中国石油化工股份有限公司 基于台站与三维地震联合采集资料的速度反演方法
CN111638551A (zh) * 2019-03-01 2020-09-08 中国石油天然气集团有限公司 地震初至波走时层析方法及装置
CN110187382A (zh) * 2019-03-05 2019-08-30 中国石油大学(华东) 一种回折波和反射波波动方程旅行时反演方法
CN112305595A (zh) * 2019-07-24 2021-02-02 中国石油化工股份有限公司 基于折射波分析地质体结构的方法及存储介质
CN112305595B (zh) * 2019-07-24 2024-05-17 中国石油化工股份有限公司 基于折射波分析地质体结构的方法及存储介质
CN112485825A (zh) * 2019-09-11 2021-03-12 中国石油化工股份有限公司 一种基于初至波走时层析的微测井解释方法
CN112485825B (zh) * 2019-09-11 2024-04-09 中国石油化工股份有限公司 一种基于初至波走时层析的微测井解释方法
CN113589374A (zh) * 2020-04-30 2021-11-02 中国石油化工股份有限公司 基于射线密度的有效速度提取方法
CN112596103A (zh) * 2020-11-24 2021-04-02 中国地质科学院地球物理地球化学勘查研究所 射线追踪方法、装置和电子设备
CN113777654A (zh) * 2021-08-06 2021-12-10 同济大学 一种基于伴随状态法初至波走时层析的海水速度建模方法
CN113791447A (zh) * 2021-10-12 2021-12-14 同济大学 一种反射结构导引的反射波层析反演方法
CN114879249A (zh) * 2022-04-13 2022-08-09 中国海洋大学 基于四面体单元走时扰动插值的地震波前走时计算方法
CN114879249B (zh) * 2022-04-13 2023-04-28 中国海洋大学 基于四面体单元走时扰动插值的地震波前走时计算方法

Similar Documents

Publication Publication Date Title
CN107843922A (zh) 一种基于地震初至波和反射波走时联合的层析成像方法
CN105277978B (zh) 一种确定近地表速度模型的方法及装置
CN102759746B (zh) 一种变偏移距垂直地震剖面数据反演各向异性参数方法
CN106094032B (zh) 一种构建地层速度模型的方法
Vinje et al. 3-D ray modeling by wavefront construction in open models
CN104216009B (zh) 一种斜井三维垂直地震剖面时间偏移的方法
CN105301636B (zh) 速度模型的建立方法和装置
CN105093319B (zh) 基于三维地震数据的地面微地震静校正方法
CN104237937B (zh) 叠前地震反演方法及其系统
CN105093292A (zh) 一种地震成像的数据处理方法和装置
CN105549077B (zh) 基于多级多尺度网格相似性系数计算的微震震源定位方法
CN105974479A (zh) Gpu空间网格体的层析2d/3d各向异性深度域速度建模方法
CN107817516A (zh) 基于初至波信息的近地表建模方法及系统
CN109001813A (zh) 一种压制多次波的方法、装置及系统
CN106338764B (zh) 生物启发计算地层圈闭油气藏超剥线识别方法
CN102338887B (zh) 不规则尺寸空变网格层析成像静校正方法
CN108303736B (zh) 各向异性ti介质最短路径射线追踪正演方法
CN106680870A (zh) 高精度地震波走时射线追踪方法
CN111638551A (zh) 地震初至波走时层析方法及装置
CN104597496B (zh) 一种二维地震资料中速度的三维空间归位方法
CN117991331B (zh) 一种基于地震监测的二维复杂模型中多次波射线追踪方法
CN105093279A (zh) 针对山前带的三维地震初至波菲涅尔体层析反演方法
CN107340537A (zh) 一种p-sv转换波叠前逆时深度偏移的方法
CN1292263C (zh) 一种用于地震勘探中射线追踪的方法
CN107942388A (zh) 一种山区地表情况下的三角网格逆时偏移方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
RJ01 Rejection of invention patent application after publication
RJ01 Rejection of invention patent application after publication

Application publication date: 20180327