CN113777654B - 一种基于伴随状态法初至波走时层析的海水速度建模方法 - Google Patents
一种基于伴随状态法初至波走时层析的海水速度建模方法 Download PDFInfo
- Publication number
- CN113777654B CN113777654B CN202110900298.1A CN202110900298A CN113777654B CN 113777654 B CN113777654 B CN 113777654B CN 202110900298 A CN202110900298 A CN 202110900298A CN 113777654 B CN113777654 B CN 113777654B
- Authority
- CN
- China
- Prior art keywords
- travel time
- sea water
- accompanying
- model
- arrival wave
- 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.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 61
- 239000013535 sea water Substances 0.000 title claims abstract description 47
- 238000004587 chromatography analysis Methods 0.000 title claims abstract description 22
- 238000004364 calculation method Methods 0.000 claims abstract description 10
- 230000005540 biological transmission Effects 0.000 claims abstract description 8
- 238000001514 detection method Methods 0.000 claims abstract description 6
- 238000003384 imaging method Methods 0.000 claims abstract description 6
- 238000007781 pre-processing Methods 0.000 claims abstract description 3
- 230000001932 seasonal effect Effects 0.000 claims abstract 2
- 238000012545 processing Methods 0.000 claims description 6
- 238000003325 tomography Methods 0.000 description 10
- 230000002159 abnormal effect Effects 0.000 description 7
- 238000005259 measurement Methods 0.000 description 5
- 238000004458 analytical method Methods 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 3
- 239000011159 matrix material Substances 0.000 description 3
- 230000007547 defect Effects 0.000 description 2
- 238000002474 experimental method Methods 0.000 description 2
- 238000000605 extraction Methods 0.000 description 2
- 238000000691 measurement method Methods 0.000 description 2
- 238000007796 conventional method Methods 0.000 description 1
- 230000003247 decreasing effect Effects 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 230000005284 excitation Effects 0.000 description 1
- 238000005286 illumination Methods 0.000 description 1
- 239000003607 modifier Substances 0.000 description 1
- 238000002945 steepest descent method Methods 0.000 description 1
Images
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/38—Seismology; Seismic or acoustic prospecting or detecting specially adapted for water-covered areas
-
- 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/30—Analysis
- G01V1/307—Analysis for determining seismic attributes, e.g. amplitude, instantaneous phase or frequency, reflection strength or polarity
-
- 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
- Life Sciences & Earth Sciences (AREA)
- Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Oceanography (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明涉及一种基于伴随状态法初至波走时层析的海水速度建模方法,包括以下步骤:1)提取OBS资料中海水透射波到达时,进行预处理得到炮点检波点和对应走时数据,即观测走时;2)根据一般情况下海水速度分布设计初始模型;3)进行基于伴随状态法初至波走时层析方法的反演迭代,得到最终成像结果。与现有技术相比,本发明具有受季节洋流因素影响小、计算效率高、适用于海洋地震勘探等优点。
Description
技术领域
本发明涉及勘探地震学中的地震成像领域,尤其是涉及一种基于伴随状态法初至波走时层析利用OBS数据的海水速度建模方法。
背景技术
初至波是地震资料中最先到达的地震信息,包括折射波和直达波,具有能量强和震相清晰的特点,能够明确地观察并拾取。目前工业界已经有利用计算机自动拾取地震资料中初至波信息的成熟方法,由于初至波信息不包括深层的反射波等信息,初至波层析成像技术目前一般应用于建立近地表速度模型,初至波信息中包括走时和波形,目前最成熟的层析成像是初至波走时层析成像方法,该方法与波动方程层析成像相比,具有更好的稳定性和更高的计算效率。
现代地震勘探中,三维地震勘探方式应用越来越广泛,地震信息数据量越来越大,基于射线追踪的传统走时层析方法依然需要巨大的计算量。谢春等2012年引入伴随状态法初至波走时层析技术,推导了基于声波方程的有限频伴随状态法层析成像,并将该方法与传统射线层析成像方法分别应用于理论模型实验和实际资料的处理,结果表明伴随状态法初至波走时层析方法可以取得与传统射线层析近似的反演结果,但计算效率得到大幅提升。2017年,李永德等提出使用近似Heissian矩阵进行预条件的伴随状态法初至波走时层析成像,该方法通过增加一次计算得到类似于射线密度的矩阵,矩阵的逆作为预条件,提高了初始模型不理想情况下的反演精度,进行理论模型和实际资料处理结果实验表明,预条件的伴随状态法初至波走时层析成像既保留了伴随状态法初至波走时层析的优点,又可以克服一阶方向的局限,提高了反演精度和效率,2021年,董良国等提出了不依赖地表法向量的改进的伴随状态法走时层析方法,克服了传统方法中伴随场依赖地表法向量的缺陷,检波点出走时残差可以正确地回传,从而提高了反演的精度。
近年来,随着海上地震勘探技术不断完善,人们对于地震资料处理精度的要求越来越高,每个参数都要求尽可能准确。目前海水中速度测量主要有两种方法:直接测量与间接测量,直接测量方法运用的仪器通常称为声速仪,一般利用收发换能器,在固定距离里测量声速;间接测量法是通过利用水文仪器测量海水的温度、盐度以及深度,然后利用环境测量与声速的经验公式,进而通过计算得到声速剖面。但以上两种常见方法在测量的过程中受到温度、盐度和压强的影响较大,并随着空间位置和时间不断变化,导致测量结果不够准确。因此,如果能用实际海上地震勘探资料来确定海水中声波速度,那么其测量精度将会大大提高,从而提高海上地震资料处理的精度和效率。
发明内容
本发明的目的就是为了克服上述现有技术存在的缺陷而提供一种基于伴随状态法初至波走时层析的海水速度建模方法。
本发明的目的可以通过以下技术方案来实现:
一种基于伴随状态法初至波走时层析的海水速度建模方法,包括以下步骤:
1)提取OBS资料中海水透射波到达时,进行预处理得到炮点检波点和对应走时数据,即观测走时;
2)根据一般情况下海水速度分布设计初始模型;
3)进行基于伴随状态法初至波走时层析方法的反演迭代,得到最终成像结果。
所述的步骤3)具体包括以下步骤:
31)使用直接扫描式最短路径法求解程函方程,计算当前模型的走时场;
32)计算当前模型的走时场与观测走时之间的走时时差,并判断目标函数是否满足终止条件,如满足则停止迭代输出反演结果,否则继续反演,进行步骤33);
33)根据走时时差和伴随方程计算伴随场λ;
34)根据伴随场λ和当前模型计算梯度;
35)确定步长,修改模型,输出当前反演结果;
36)重复步骤31-35),直到获得最终的反演结果。
所述的步骤31)具体包括以下步骤:
311)确定震源位置并确定震源相邻节点的最小走时;
312)直接在各方向上进行全场扫描,得到当前模型的走时场。
所述的步骤33)具体包括以下步骤:
331)根据当前模型的走时场与实际数据提取到的海水透射波走时做差值;
332)根据伴随方程计算得到检波点处的伴随变量;
333)根据伴随方程扫描得到全空间的伴随场。
所述的步骤332)中,伴随方程的表达式:
其中,λ(x)为当前位置x处的伴随变量,x为当前位置,Ω为计算区域,t(x)为当前模型的走时场,Tobs(x)为观测得到的走时场,δ为δ函数,且表示检波点序号,ri为第i个检波点位置,N为检波点数量,/>为梯度算子。
所述的步骤34)中,伴随变量和梯度满足关系式:
其中,J(x)为目标函数,c为模型速度。
所述的目标函数的表达式为:
所述的步骤2)中,在海洋地震数据处理过程中,一般情况下假设海水速度为常数,当测区存在更详细的海水速度分布信息时,则在此基础上优化初始模型。
假设海水速度为1500米/秒。
所述的步骤1)中,OBS资料来源于在短时间内一次海上主动地震采集到的OBS走时数据,以此降低季节和洋流因素的影响。
与现有技术相比,本发明具有以下优点:
一、本发明根据在短时间内一次海上主动地震采集到的大量OBS走时数据,得到的走时数据,受到的季节、洋流等因素影响较小。
二、由于引入了伴随状态法,本发明中的梯度求取过程避免了Fréchet导数的求取,而是间接计算得到梯度,这样的计算效率更高,易于实现并行计算。
三、本发明中的海水速度反演流程与地震勘探反演流程相似,在实际海洋地震勘探过程容易实现。
附图说明
图1为本发明的方法流程图。
图2为实施例1中的理论海水速度模型。
图3为实施例1中的反演初始模型。
图4为实施例1中的透射波走时层析成像结果。
图5为实施例1中的反演结果竖直方向速度抽线,其中,图(5a)为X=500米处的速度抽线,图(5b)为X=750米处的速度抽线。
图6为实施例1中的反演结果水平方向速度抽线,其中,图(6a)为Y=200米处的速度抽线,图(6b)为Y=700米处的速度抽线。
图7为实施例2中的理论海水速度模型。
图8为实施例2中的反演初始模型。
图9为实施例2中的透射波走时层析成像结果。
图10为实施例2中的反演结果速度抽线,其中,图(10a)为X=600米,Y=600米处的垂直方向速度抽线,图(10b)为Z=250米,Y=600米处的水平方向速度抽线。
图11为实施例3中的反演初始模型。
图12为实施例3中的透射波走时层析成像结果。
具体实施方式
下面结合附图和具体实施例对本发明进行详细说明。
实施例1:
根据海水速度分布特性,本例中设置以下二维理论模型:模型大小为2000米*1000米,网格为200*100,间距为10米,背景速度变化趋势为“高速-低速-高速”,最大速度为1500米/秒,在500m左右的深度时速度达到最小,在(1000,600)位置设置一个高速异常体,速度为1480米/秒,950米以下为海底,速度为3000米/秒,如图2所示。将炮点设置到模型顶部即海面,OBS设置于海底,炮点个数201炮,炮间距为10米,OBS个数为201个,间距为10米。使用海水速度为1500m/s的常速度模型作为初始模型进行反演,初始模型如图3所示。初始步长为50,随目标函数下降。反演迭代603代,初始目标函数为39.9,最终目标函数为0.0066,目标函数下降到初始模型情况下的0.016%,总共用时1小时31分。
具体实施方式如下:
1)通过直接扫描式最小路径法对理论模型进行地震波射线传播正演,得到理论走时;
2)通过直接扫描式最小路径法对初始模型进行地震波射线传播正演,得到计算走时;
3)计算接收点处的理论走时和计算走时的差,通过伴随方程计算得到伴随变量场;
4)根据伴随变量与梯度的关系计算得到梯度场;
5)根据最速下降法,确定模型修改量,得到新速度模型;
6)将新速度模型作为初始模型继续迭代,重复步骤2-5),直到新速度模型满足收敛条件,即可结束迭代,输出最终结果
得到反演结果如图4所示。对反演结果进行观察,可以看出改进的伴随状态法走时层析可以比较准确地反演出海水的速度分布情况,对于高速异常体位置和大小有比较好的表达,走时基本拟合,边界处误差较为明显。水平方向分辨率明显高于垂直方向,垂直方向上异常体的形状发生了拉伸,背景速度趋于平均。
分别在X=500米,X=750米,Y=200米,Y=700米处对反演结果分别进行水平方向上和垂直方向上速度抽线分析,可以很明显的看出水平方向上的反演精度高,分辨率高,既能很好地趋近真实速度,也能很好地表达速度横向变化;垂向上反演分辨率较低,反演精度较低,反演速度可以大致趋向真实速度,出现水平方向和垂直方向分辨率差别很大的原因是由于采用了单边激发单边接收的观测系统和走时反演,这样的条件下,反演对于沿射线方向的速度变化不敏感,对于射线之间的速度变化敏感。
实施例2:
根据以上二维理论模型,设计了三维理论模型,模型大小为1000米*1000米*500米,网格为100*100*50,网格间距10米,海水速度部分设计沿用二维的思路,海水背景速度呈“高-低-高”变化的趋势,最大速度1500米/秒,到200米深度速度达到最小值,模型中心位置存在一个高速异常体,异常体速度为1480米/秒,如图6所示。为模拟OBS观测系统,将炮点设置于模型顶部即海面,炮数为51*51共2601炮,炮间距10米,均匀布满海面;OBS位于模型底部即海底,个数为101*101共10201个,间距10米,均匀布满海底。使用海水速度为1500米/秒的均匀模型作为初始模型进行反演,初始模型如图7所示。反演迭代21轮,初始目标函数为35.5,最终目标函数为0.87,目标函数下降到初始模型的2.4%,具体实施流程与实例1相似,区别在于将二维情况的计算拓展至三维。
对反演结果进行观察,可以看出三维情况下改进的伴随状态法走时层析可以比较准确地反演出海水的速度分布情况,对于高速异常体位置和大小有比较好的表达,走时基本拟合,水平方向分辨率明显高于垂直方向,垂直方向上异常体的形状发生了拉伸,背景速度趋于平均。
在X=600米,Y=600米处进行速度垂向抽线分析,在Y=600米,Z=250米处进行速度水平抽线分析,可以很明显的看出水平方向上的反演精度高,分辨率高,既能很好地趋近真实速度,也能很好地表达速度横向变化;垂向上反演分辨率较低,反演精度较低,反演速度可以大致趋向真实速度。
实施例3:
本实施例将本发明提出的基于OBS的海水速度建模方法应用于我国东海的真实OBS数据,反演使用的初始模型如图11所示,初始模型水平长度11000米,深度90米;网格大小为11000*9,间距1米,海水速度为均匀1500米/秒。共91个OBS,深度为81米,分布从638米到10715米;炮点深度9米,炮点分布从0米到10999米,共283炮。海底平均深度85米。反演迭代14次后停止,初始目标函数为788.2,最终目标函数为226.3,目标函数下降到使用初始模型时的28.7%。
分析反演结果,可以看出改进的伴随状态法初至波走时层析成像方法应用与实际资料时,可以大致表现出真实海水的速度分布,反演结果符合海水理论速度分布;同时水平方向上对于海水速度的变化有较好表现,垂直方向上海水速度变化较小,推测原因其一是因为实际资料测区海水较浅,实际海水垂向上变化较小;其二是因为观测系统的局限性,使得垂向分辨率低于横向分辨率。此外,计算区域两侧明显反演精度比较低,推测原因是该区域射线密度比较稀疏,照明不足。
Claims (6)
1.一种基于伴随状态法初至波走时层析的海水速度建模方法,其特征在于,包括以下步骤:
1)提取OBS资料中海水透射波到达时,进行预处理得到炮点检波点和对应走时数据,即观测走时;
2)根据一般情况下海水速度分布设计初始模型;
3)进行基于伴随状态法初至波走时层析方法的反演迭代,得到最终成像结果;
所述的步骤3)具体包括以下步骤:
31)使用直接扫描式最短路径法求解程函方程,计算当前模型的走时场;
32)计算当前模型的走时场与观测走时之间的走时时差,并判断目标函数是否满足终止条件,如满足则停止迭代输出反演结果,否则继续反演,进行步骤33);
33)根据走时时差和伴随方程计算伴随场λ;
34)根据伴随场λ和当前模型计算梯度;
35)确定步长,修改模型,输出当前反演结果;
36)重复步骤31-35),直到获得最终的反演结果;
所述的步骤31)具体包括以下步骤:
311)确定震源位置并确定震源相邻节点的最小走时;
312)直接在各方向上进行全场扫描,得到当前模型的走时场;
所述的步骤33)具体包括以下步骤:
331)根据当前模型的走时场与实际数据提取到的海水透射波走时做差值;
332)根据伴随方程计算得到检波点处的伴随变量;
333)根据伴随方程扫描得到全空间的伴随场;
所述的步骤332)中,伴随方程的表达式:
4.根据权利要求1所述的一种基于伴随状态法初至波走时层析的海水速度建模方法,其特征在于,所述的步骤2)中,在海洋地震数据处理过程中,一般情况下假设海水速度为常数,当测区存在更详细的海水速度分布信息时,则在此基础上优化初始模型。
5.根据权利要求4所述的一种基于伴随状态法初至波走时层析的海水速度建模方法,其特征在于,假设海水速度为1500米/秒。
6.根据权利要求1所述的一种基于伴随状态法初至波走时层析的海水速度建模方法,其特征在于,所述的步骤1)中,OBS资料来源于在短时间内一次海上主动地震采集到的OBS走时数据,以此降低季节和洋流因素的影响。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202110900298.1A CN113777654B (zh) | 2021-08-06 | 2021-08-06 | 一种基于伴随状态法初至波走时层析的海水速度建模方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202110900298.1A CN113777654B (zh) | 2021-08-06 | 2021-08-06 | 一种基于伴随状态法初至波走时层析的海水速度建模方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN113777654A CN113777654A (zh) | 2021-12-10 |
CN113777654B true CN113777654B (zh) | 2023-07-04 |
Family
ID=78836858
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202110900298.1A Active CN113777654B (zh) | 2021-08-06 | 2021-08-06 | 一种基于伴随状态法初至波走时层析的海水速度建模方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN113777654B (zh) |
Family Cites Families (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US10067255B2 (en) * | 2015-09-04 | 2018-09-04 | Saudi Arabian Oil Company | Automatic quality control of seismic travel time |
CN109541681B (zh) * | 2017-09-22 | 2020-07-17 | 中国海洋大学 | 一种拖缆地震数据和少量obs数据联合的波形反演方法 |
WO2019071504A1 (zh) * | 2017-10-12 | 2019-04-18 | 南方科技大学 | 一种基于两点射线追踪的地震走时层析反演方法 |
CN107843922A (zh) * | 2017-12-25 | 2018-03-27 | 中国海洋大学 | 一种基于地震初至波和反射波走时联合的层析成像方法 |
CN109100798B (zh) * | 2018-06-22 | 2019-10-25 | 广州海洋地质调查局 | 实现折射多次波层析反演的方法、装置及处理终端 |
CN109085643A (zh) * | 2018-07-30 | 2018-12-25 | 中国石油化工股份有限公司 | 早至波的分步联合反演方法 |
CN111781639B (zh) * | 2020-06-04 | 2021-06-04 | 同济大学 | 针对obs多分量数据的炮检互易弹性波全波形反演方法 |
CN112363211A (zh) * | 2020-11-23 | 2021-02-12 | 同济大学 | 一种改进的sirt法射线走时层析成像方法 |
-
2021
- 2021-08-06 CN CN202110900298.1A patent/CN113777654B/zh active Active
Also Published As
Publication number | Publication date |
---|---|
CN113777654A (zh) | 2021-12-10 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN102027390B (zh) | 对地震数据的联合内插及消除虚反射 | |
US7068568B2 (en) | Method of processing seismic data | |
CN106526674B (zh) | 一种三维全波形反演能量加权梯度预处理方法 | |
CN106950568B (zh) | 一种自适应多节点等效声速剖面的构建方法 | |
CN106405651B (zh) | 一种基于测井匹配的全波形反演初始速度模型构建方法 | |
MX2011006036A (es) | Uso de inversion de forma de onda para determinar las propiedades de un medio en el subsuelo. | |
AU2001274411A1 (en) | A method of processing seismic data | |
MX2014014770A (es) | Analisis de datos sismicos utilizando coleccion de datos de nodos del fondo oceanico. | |
CN113552625B (zh) | 一种用于常规陆域地震数据的多尺度全波形反演方法 | |
WO2008067588A1 (en) | Method of building a sub surface velocity model | |
CN111487678B (zh) | 一种确定高分辨率小多道地震最小偏移距和系统延迟的分析方法 | |
CN104374828A (zh) | 一种隐患探测的超声波层析成像方法 | |
CN107229071B (zh) | 一种地下构造反演成像方法 | |
CN102901985A (zh) | 一种适用于起伏地表的深度域层速度修正方法 | |
CN109633749A (zh) | 基于散射积分法的非线性菲涅尔体地震走时层析成像方法 | |
WO2018145487A1 (zh) | 一种测量海水中声波速度的方法 | |
CN112099082B (zh) | 一种共面元共方位角道集的地震回折波走时反演方法 | |
CA2412995C (en) | Seismic survey system | |
Gong et al. | Combined migration velocity model-building and its application in tunnel seismic prediction | |
CN109782355A (zh) | Obs检波点漂移的检测方法及装置 | |
CN113777654B (zh) | 一种基于伴随状态法初至波走时层析的海水速度建模方法 | |
CN109490961B (zh) | 起伏地表无射线追踪回折波层析成像方法 | |
CN113568041B (zh) | 时移地震三维拖缆采集数据的可重复性分析方法及系统 | |
CN108375794B (zh) | 基于对称观测的vsp缝洞绕射成像技术方法 | |
CN112257241B (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 |