CN109541691A - 一种地震速度反演方法 - Google Patents

一种地震速度反演方法 Download PDF

Info

Publication number
CN109541691A
CN109541691A CN201811460593.4A CN201811460593A CN109541691A CN 109541691 A CN109541691 A CN 109541691A CN 201811460593 A CN201811460593 A CN 201811460593A CN 109541691 A CN109541691 A CN 109541691A
Authority
CN
China
Prior art keywords
big gun
super big
super
wavelet
velocity field
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.)
Granted
Application number
CN201811460593.4A
Other languages
English (en)
Other versions
CN109541691B (zh
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.)
China University of Petroleum East China
Original Assignee
China University of Petroleum East 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 China University of Petroleum East China filed Critical China University of Petroleum East China
Priority to CN201811460593.4A priority Critical patent/CN109541691B/zh
Publication of CN109541691A publication Critical patent/CN109541691A/zh
Application granted granted Critical
Publication of CN109541691B publication Critical patent/CN109541691B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

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/30—Analysis
    • G01V1/303—Analysis for determining velocity profiles or travel times
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/62—Physical property of subsurface
    • G01V2210/622—Velocity, density or impedance
    • G01V2210/6222—Velocity; travel time

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

本发明公开了一种地震速度反演方法，包括：输入初始的观测地震数据、初始的速度模型与极性编码矩阵；对观测炮集利用极性编码矩阵合成编码超级炮集；利用编码超级炮集与极性编码矩阵生成加权编码的参考道；利用参考道和目标子波对超级炮记录采用维纳滤波生成目标子波对应的超级炮；利用目标子波和编码函数与速度场生成正演超级炮；利用正演超级炮与生成的超级炮记录做差反传，得到反演梯度；利用抛物插值计算最佳步长，并对速度场进行更新；重复S2到S7直到满足迭代次数或者速度场要求，输出最终的速度场。本发明对维纳滤波采用的参考道是一种加权的多震源地震道，其不需要估计子波，减少子波估计的计算量。

Description

一种地震速度反演方法
技术领域
本发明涉及一种地震速度反演方法，属于油气物探工程领域。
背景技术
在当前的勘探领域中，地震速度场的反演对于地震数据的偏移成像和解释具有重要的意义，全波形反演能够充分利用地震数据中所包含的振幅、相位、走时、波形等信息，对地下模型参数进行精细刻画。
在当前，通过维纳滤波实现多尺度的全波形反演是一种常用的手段。但是在这种方式下，维纳滤波必须使用原始子波，存在原始子波依赖性强、计算量巨大、反演稳定性等诸多问题，由于采集成本、坏道、噪声、地形等诸多因素的影响，实际采集到的数据往往提取的原始子波不准确，给地震数据的速度反演造成了一定的困扰。
因此，发展一种不依赖子波的快速地震速度反演方法是十分必要的。
发明内容
本发明的目的在于，提供一种地震速度反演方法，它可以解决当前技术中存在的问题，对不知道地震子波的地震数据，提高地震数据反演速度场的精度。
为解决上述技术问题，本发明提供一种地震速度反演方法，包括以下步骤：
S1，输入初始的观测地震数据、初始的速度模型与极性编码矩阵；
S2，对观测炮集利用极性编码矩阵合成编码超级炮集；
S3，利用编码超级炮集与极性编码矩阵生成加权编码的参考道；
S4，利用参考道和目标子波对超级炮记录采用维纳滤波生成目标子波对应的超级炮；
S5，利用目标子波和编码函数与速度场生成正演超级炮；
S6，利用正演超级炮与S4生成的超级炮记录做差反传，互相关得到反演梯度；
S7，利用抛物插值计算最佳步长，并对速度场进行更新；
S8，重复S2到S7直到满足迭代次数或者速度场要求，输出最终的速度场。
与现有技术相比，本发明对维纳滤波采用的参考道是一种加权的多震源地震道，其不需要估计子波，减少子波估计的计算量。
附图说明
附图是用来提供对本发明的进一步理解，并且构成说明书的一部分，与下面的具体实施方式一起用于解释本发明，但并不构成对本发明的限制。在附图中：
图1为本发明实施例的流程示意图；
图2为本发明的编码矩阵示意图；
图3为本发明使用的速度场模型图；
图4为本发明使用的单炮记录；
图5为本发明实施例的初始速度场示意图；
图6为本发明的观测数据中的合成超级炮记录示意图；
图7为使用本发明方法生成的参考道；
图8为使用本发明得到目标子波对应的超级炮记录；
图9为使用本发明得到反演速度场；
图10为不准确子波得到反演速度场；
图11为常规不依赖子波方法得到反演速度场；
图12为准确子波方法得到反演速度场。
下面结合附图和具体实施方式对本发明作进一步的说明。
具体实施方式
本发明的实施例1：一种地震速度反演方法，如图1所示，包括以下步骤：
S1，输入初始的观测地震数据、初始的速度模型与极性编码矩阵。在实践中，根据需要，观测数据为单震源地震数据，极性编码矩阵如图2所示，其中编码矩阵是一个K×Ns的矩阵，其中值为随机生成的正负1，K为总迭代次数，Ns为编码的总炮数，由于每次迭代采用的编码极性是随机的，因此其可以改善不同炮记录间的差异性，较好的压制串扰噪音；为了说明本发明方法的有效性例如，将本发明的方法应用到国际标准Marmous模型上如图3所示，对图3模型，使用有限差分正演得到很多数量的单炮记录，将其作为观测数据，选取炮点为中间位置的单炮数据展示，如图4；图5为所述发明方法利用观测记录进行反演使用初始的速度模型。
S2，对观测炮集利用极性编码矩阵合成编码超级炮集。本文反演速度方法选取的是全波形反演方法，其是一种迭代反演方法，对S1生成的炮集(如图4)利用极性编码矩阵(图2)进行编码形成的超级炮集，则第k次迭代合成的超级炮记录Usu,k，具体公式如下：
其中，uobs,n表示炮集第n炮的观测炮记录，en,k表示极性编码矩阵第k次迭代第n炮的极性值，其中编码矩阵是一个K×Ns的矩阵(如图2)，其中值为随机生成的正负1，K为总迭代次数，Ns为编码的总炮数，Usu,k为第k次迭代合成的超级炮，因为其将很多炮记录合成了一个超级炮，因此计算量大大减少。因为采用极性编码方式，加大了临炮之间的相干性，其可以较好的压制多震源间串扰噪音，更精确的对复杂区域进行速度场的反演。对S2的所有单炮数据(图4)利用极性编码的公式进行编码，其中得到第k次迭代的超级炮记录Usu,k，如图6所示。
编码超级炮集中包含多个超级炮记录。在进行迭代时，可以先根据部分或者全体观测炮记录生成多个的超级炮，组成超级炮集。在这个过程中，观测炮集不变，en,k随着迭代次数进行将会不同。本发明中，每次迭代，可以使用所有单炮记录生成了一个超级炮记录。
S3，利用编码超级炮集与极性编码矩阵生成加权编码的参考道。利用S2中合成的编码超级炮集Usu,k(图6)与S1中所用的极性编码矩阵(图2)生成加权编码的参考道，具体公式如下：
Usu,k(x(n))为S2中第k次迭代合成的超级炮，x(n)第n个炮记录对应的炮点位置，Uref,k为第k次迭代生成所需要的参考道，如图7展示了本发明方法第1次迭代生成参考道。由于参考道选取是对合成的超级炮进行处理，其相对于原始所有的单炮集(图4)，存储较小，参考道选取也是相乘相加运算，计算简单，耗时较小。由于采用了编码矩阵相乘，相当于对超级炮是一个加权处理，通过加权处理其可以与原始子波更加接近，每次迭代加权系数不一样，这样对于不同于原始子波的部分，可以在迭代中消除误差的影响，参考道选取也是相乘相加运算，计算简单，耗时较小。
S4，利用参考道和目标子波对超级炮记录采用维纳滤波生成目标子波对应的超级炮。
在多尺度全波形反演方法利用维纳滤波是一种比较常见的方法。在利用维纳滤波进行多震源多尺度反演的时候，需要给定一个原始的地震子波和主频为w的目标子波，其中目标子波一般可以选取主频为w的雷克子波，但是原始子波的振幅、相位等信息一般不知道，需要反演或者计算得到，如果求取的原始子波不准确，则影响反演结果的精度。应用维纳滤波公式，本文的主频为w的雷克子波对应的超级炮生成具体公式：
Star(w)为主频为w的雷克子波，Utar,k(w)为第k次迭代目标子波Star(w)所对应的超级炮，Uref,k为应用S3得到第k次迭代的参考道，Usu,k为应用S2公式生成的第k次迭代超级炮，ε为远小于1的正数。利用S3所求取的参考道(图7)Uref,k和S2生成的超级炮Usu,k，以及主频为w的雷克子波，利用S4公式生成目标子波Star(w)对应的第k次迭代的超级炮集，如图8所示。
应用维纳滤波公式，本发明中目标子波为主频为w的雷克子波，原始子波用S3生成的参考道Uref，k代替，避免了求取原始子波的过程，减少子波估计的计算量。
S5，利用目标子波和编码矩阵与速度场生成正演超级炮。利用S1中的极性编码矩阵(图2)和主频为w的雷克子波Star(w)生成超级炮震源，对第k-1次迭代得到速度场模型vk-1利用有限差分进行正演生成第k次的超级炮Dcal,k，如果k＝1，则vk-1为S1中的初始速度场模型(图5)具体公式如下：
L为有限差分正演算子，v为初始的速度场，xn为第n炮的炮点位置，en,k表示极性编码矩阵(如图2)第k次迭代第n炮的极性值，其为随机生成的正负1。
S6，利用S5中正演超级炮Dcal,k与S4生成的超级炮记录Utar,k(w)(图8)做差，并进行波场反传，与正传波场互相关得到反演梯度。反演速度场的过程即为通过不断迭代最小化一个目标函数，这里的目标函数即正演超级炮与S4生成的超级炮记录，表示为：
χ(v,w)＝||Dcal,k(w)-Utar,k(w)||2 (1)
其中，χ(v)表示速度场v对应的数据残差二范数。根据数据残差确定梯度，基于梯度进行速度场模型迭代，具体的说，对数据残差进行反传得到第k次的更新梯度值gk，梯度更新公式如下：
其中，vk表示第k次迭代的反演速度场，δ为求导算子，w为S4第k次迭代所用的雷克子波的主频。由于本发明方法采用的超级炮反演方法，其计算效率相比普通单震源反演大大提高。
S7，利用抛物插值计算最佳步长，并对速度场进行更新。用线性搜索方法或者抛物拟合方法等方法可以求取步长，使得第k次迭代反演速度场对应的S6中的误差最小，本文利用抛物插值方法求取第k次迭代的更新步长αk，然后利用S6中的梯度更新速度场:
vk＝vk-1-αkgk
其中，vk表示第k次迭代的反演速度场，gk为S6中求取的梯度场，vk-1为第k-1次迭代的反演速度场，如果是第一次迭代k＝1，则vk-1为S1初始速度场模型(图5)。
S8，重复S2到S7直到满足迭代次数或者速度场要求，输出最终的速度场。当误差小于预设值时，所得的速度场模型即为最终的速度场模型。在对S1中Marmousi模型(图3)的反演中，使用初始速度场模型(图5)经过反复迭代，应用本文的发明方法得到第120次迭代的反演结果如图9所示，可以看出图9中对背景场及复杂构造层位反演更好，可以较为准确的反演速度场，减少子波不准确的影响。如果采用不准确子波作为原始子波，以图5作为初始的反演速度，以S2中超级炮作为输入炮集，采用维纳滤波的多尺度全波形反演方法，第120次迭代的多震源反演速度场结果如图10，可以看出反演的结果不准确。如果使用第1炮位置下面的地震道作为参考道，以图5作为初始的反演速度，以S2中超级炮作为输入炮集，进行维纳滤波的多震源全波形反演，第120次迭代的反演结果如图11，可以看出反演的速度场相对于本文的反演方法不太准确。如果使用已知的准确的地震子波，S1中有限差分正演得到观测记录(图4)的地震子波，以图5作为初始的反演速度，以S2中超级炮作为输入炮集，进行维纳滤波的多震源全波形反演，第120次迭代的反演速度场结果如图12，可以看出反演的速度场相对于本文的反演方法基本一致。但是对于实际的观测记录我们比较难得到准确原始子波，选取不同方法提取的子波时得到的原始子波也不一样，不一定匹配对应的地震记录，而本文的方法不用求取原始子波，参考道选取也是相乘相加运算，计算简单，耗时较小。
通过前述的方法，发明对维纳滤波采用的参考道是一种加权的多震源地震道，其不需要估计子波，减少子波估计的计算量，提高了效率。此外，本发明对观测数据和震源采用极性编码矩阵，充分加大混叠炮之间的差异性，减少临炮之间的干扰，在混叠炮数据反演时较为好的压制串扰噪音，提高了多震源数据速度场反演时的精度。

Claims (4)

1.一种地震速度反演方法，其特征在于，包括以下步骤：
S1，输入初始的观测地震数据、初始的速度模型与极性编码矩阵；
S2，对观测炮集利用极性编码矩阵合成编码超级炮集，其中，编码超级炮集包含多个超级炮记录；
S3，利用编码超级炮集与极性编码矩阵生成加权编码的参考道；
S4，利用参考道和目标子波对超级炮记录采用维纳滤波生成目标子波对应的超级炮；
S5，利用目标子波和编码矩阵与速度场生成正演超级炮；
S6，利用正演超级炮与S4生成的超级炮记录做差反传，互相关得到反演梯度；
S7，利用抛物插值计算最佳步长，并对速度场进行更新；
S8，重复S2到S7直到满足迭代次数或者速度场要求，输出最终的速度场。
2.根据权利要求1所述的方法，其特征在于，所述步骤S2合成编码超级炮集，具体公式如下：
其中，uobs,n表示炮集第n炮的炮记录，en,k表示极性编码矩阵第k次迭代第n炮的极性值，其中编码矩阵是一个K×Ns的矩阵，其中值为随机生成的正负1,K为总迭代次数，Ns为编码的总炮数，Usu,k为第k次迭代合成的超级炮。
3.根据权利要求2所述的方法，其特征在于，所述步骤S3中的编码超级炮集与极性编码矩阵生成加权编码的参考道，具体公式如下：
Usu,k(x(n))为S2中第k次迭代合成的超级炮，x(n)第n个炮记录对应的炮点位置，Uref,k为第k次迭代生成的加权编码的参考道。
4.根据权利要求3所述的一种地震速度反演方法，其特征在于，所述步骤S4中，利用参考道和目标子波对超级炮记录采用维纳滤波生成目标子波对应的超级炮，具体公式：
Star(w)为主频为w的目标子波，Utar,k(w)为第k次迭代目标子波Star(w)所对应的超级炮，Uref,k为应用S3得到第k次迭代的参考道，Usu,k为应用S2公式生成的第k次迭代超级炮，ε为小于1的正数。
CN201811460593.4A 2018-12-01 2018-12-01 一种地震速度反演方法 Expired - Fee Related CN109541691B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201811460593.4A CN109541691B (zh) 2018-12-01 2018-12-01 一种地震速度反演方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201811460593.4A CN109541691B (zh) 2018-12-01 2018-12-01 一种地震速度反演方法

Publications (2)

Publication Number Publication Date
CN109541691A true CN109541691A (zh) 2019-03-29
CN109541691B CN109541691B (zh) 2020-05-19

Family

ID=65851742

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201811460593.4A Expired - Fee Related CN109541691B (zh) 2018-12-01 2018-12-01 一种地震速度反演方法

Country Status (1)

Country Link
CN (1) CN109541691B (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110376642A (zh) * 2019-06-02 2019-10-25 中国石油大学(华东) 一种基于锥面波的三维地震速度反演方法
CN113050179A (zh) * 2021-03-11 2021-06-29 中国科学院地质与地球物理研究所 一种三维多源探地雷达设备及方法
CN114063160A (zh) * 2020-08-10 2022-02-18 中国石油化工股份有限公司 一种地震速度反演方法及装置

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105319581A (zh) * 2014-07-31 2016-02-10 中国石油化工股份有限公司 一种高效的时间域全波形反演方法
CN107229071A (zh) * 2017-05-25 2017-10-03 中国石油大学(华东) 一种地下构造反演成像方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105319581A (zh) * 2014-07-31 2016-02-10 中国石油化工股份有限公司 一种高效的时间域全波形反演方法
CN107229071A (zh) * 2017-05-25 2017-10-03 中国石油大学(华东) 一种地下构造反演成像方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
CHAIWOOT BOONYASIRIWAT ET AL.: "An efficient multiscale method for time-domain waveform tomography", 《GEOPHYSICS》 *
YUNDONG GUO ET AL.: "Fast building initial velocity modeling using encoding multiscale multi-shot full-waveform inversion", 《SEG INTERNATIONAL EXPOSITION AND 87TH ANNUAL MEETING》 *
国运东等: "基于极性振幅编码的全波形反演策略分析", 《中国地球科学联合学术年会2015》 *

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110376642A (zh) * 2019-06-02 2019-10-25 中国石油大学(华东) 一种基于锥面波的三维地震速度反演方法
CN114063160A (zh) * 2020-08-10 2022-02-18 中国石油化工股份有限公司 一种地震速度反演方法及装置
CN114063160B (zh) * 2020-08-10 2023-03-31 中国石油化工股份有限公司 一种地震速度反演方法及装置
CN113050179A (zh) * 2021-03-11 2021-06-29 中国科学院地质与地球物理研究所 一种三维多源探地雷达设备及方法

Also Published As

Publication number Publication date
CN109541691B (zh) 2020-05-19

Similar Documents

Publication Publication Date Title
CN108873066B (zh) 弹性介质波动方程反射波旅行时反演方法
CN109407151B (zh) 基于波场局部相关时移的时间域全波形反演方法
CN106526674B (zh) 一种三维全波形反演能量加权梯度预处理方法
AU2011248989B2 (en) Artifact reduction in method of iterative inversion of geophysical data
CN107843925B (zh) 一种基于修正相位的反射波波形反演方法
CN109541691A (zh) 一种地震速度反演方法
CN105467442B (zh) 全局优化的时变稀疏反褶积方法及装置
CN110187384B (zh) 贝叶斯时移地震差异反演方法及装置
CN111239806B (zh) 基于振幅增量编码的时间域全波形反演方法
CN105093278A (zh) 基于激发主能量优化算法的全波形反演梯度算子的提取方法
WO2016005597A1 (en) Method for obtaining estimates of a model parameter so as to characterise the evolution of a subsurface volume
CN111580163A (zh) 一种基于非单调搜索技术的全波形反演方法及系统
CN111999764A (zh) 基于时频域目标函数的盐下构造最小二乘逆时偏移方法
CN110737018A (zh) Vsp地震资料各向异性建模方法
CN111208568B (zh) 一种时间域多尺度全波形反演方法及系统
CN113791447B (zh) 一种反射结构导引的反射波层析反演方法
CN110376642B (zh) 一种基于锥面波的三维地震速度反演方法
CN112305615B (zh) 一种地震资料角度域共成像点道集提取方法及系统
CN114063160B (zh) 一种地震速度反演方法及装置
CN111190224A (zh) 一种基于三维地震波反向照明的动态采样全波形反演系统及方法
CN111722287B (zh) 一种基于渐进数据同化方法的震相特征识别波形反演方法
Zhou et al. Anisotropic model building with well control
CN114185090B (zh) 岩相和弹性参数同步反演方法、装置、电子设备及介质
CN113050162B (zh) 基于Zoeppritz方程的粘弹介质地震反演方法
CN115993648A (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
GR01 Patent grant
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20200519

Termination date: 20201201

CF01 Termination of patent right due to non-payment of annual fee