CN115951407B - 多次波成像角度域共成像点道集计算方法及计算设备 - Google Patents

多次波成像角度域共成像点道集计算方法及计算设备 Download PDF

Info

Publication number
CN115951407B
CN115951407B CN202211121320.3A CN202211121320A CN115951407B CN 115951407 B CN115951407 B CN 115951407B CN 202211121320 A CN202211121320 A CN 202211121320A CN 115951407 B CN115951407 B CN 115951407B
Authority
CN
China
Prior art keywords
field
wave field
wave
source
wavefield
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
Application number
CN202211121320.3A
Other languages
English (en)
Other versions
CN115951407A (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.)
Sun Yat Sen University
Original Assignee
Sun Yat Sen University
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 Sun Yat Sen University filed Critical Sun Yat Sen University
Priority to CN202211121320.3A priority Critical patent/CN115951407B/zh
Publication of CN115951407A publication Critical patent/CN115951407A/zh
Application granted granted Critical
Publication of CN115951407B publication Critical patent/CN115951407B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/30—Assessment of water resources

Abstract

本发明涉及地震波油气勘探领域，公开了多次波成像角度域共成像点道集的计算方法及计算设备，包括：正向延拓获得多次波震源波场的位移场；根据能流密度的计算公式获得震源波场的能流密度矢量；扫描每一时刻震源波场能流密度矢量的模值，保存模值最大的震源波场能流密度矢量；逆向延拓获得多次波检波波场的位移场；根据能流密度的计算公式获得多次波检波波场的能流密度矢量；扫描每一时刻获得检波波场能流密度矢量的模值，保存模值最大的多次波检波波场能流密度矢量；计算地震波场的入射角；计算多次波成像角度域共成像点道集；本发明利用能流密度矢量计算多次波成像的角度域共成像点道集，计算准确度高，具有明确的物理意义。

Description

多次波成像角度域共成像点道集计算方法及计算设备
技术领域
本发明涉及地震波油气勘探领域，特别是涉及多次波成像角度域共成像点道集的计算方法。
背景技术
在地震勘探中，叠前深度偏移成像技术利用地下传播的地震波获得地下结构信息，角度域共成像点道集(简称角道集)作为叠前深度偏移成像技术的产物，是地震勘探中用于构建速度场、属性分析以及优化成像结果得重要工具。常规角度域共成像点道集是由叠前深度偏移技术利用一次反射的地震波(一次波)进行计算的，而对于海洋勘探，尤其是以拖缆形式的勘探采集得到的地震数据含有丰富的多次反射的地震波(多次波)，通常这种多次波在地震数据处理中往往作为噪音进行去除。
作为一次波的补充，多次波可以被用来成像而不是作为噪声。相比于一次波成像而言，多次波成像可以提供额外的角度照明，从而提高地下成像结果的分辨率，目前已在工业界得以广泛的应用(Lu,S.,2021,Migration using sea surface-related multiples:Challenges and opportunities:Geophysics,86(5),1-42)。对于多次波成像的角道集计算而言，目前主要的方法还是利用地下偏移距道集转角道集的方法(Sava,P.,和GuittonA.,2005,Multiple attenuation in the image space:Geophysics,70(195).V10–V20)。但是本申请的发明人发现：利用地下偏移距道集转角道集的方法计算多次波的角度域共成像点道集不具有明确的物理意义，因此计算的多次波角度域共成像点道集不准确。
地震波场的能流密度矢量可直接反映地震波场的传播方法，具有明确的物理意义，因此本专利利用能流密度矢量计算多次波成像的角道集。
发明内容
本发明的目的在于提供多次波成像角度域共成像点道集的计算方法及计算设备，利用能流密度矢量计算多次波成像的角度域共成像点道集，具有明确的物理意义，且保证多次波成像角度域共成像点道集的计算准确度高。
为实现上述目的，本发明的一个方面提供了多次波成像角度域共成像点道集的计算方法，其包括以下步骤：
步骤一，基于声波波动方程，输入地表记录的炮集数据作为震源，进行地震波场的正向延拓获得多次波震源波场的位移场；
步骤二，将步骤一保存的多次波震源波场位移场进行梯度运算和关于传播时间的导数运算，并根据能流密度的计算公式获得震源波场的能流密度矢量；
步骤三，在给定的时间窗内，扫描每一时刻步骤二获得的震源波场能流密度矢量的模值，并保存模值最大的震源波场能流密度矢量；
步骤四，基于声波波动方程，输入地表记录的炮集数据作为边值条件，进行地震波场的逆向延拓获得多次波检波波场的位移场；
步骤五，将步骤四保存的多次波检波波场位移场进行梯度运算和关于传播时间的导数运算，并根据能流密度的计算公式获得多次波检波波场的能流密度矢量；
步骤六，在给定的时间窗内，扫描每一时刻步骤五获得的检波波场能流密度矢量的模值，并保存模值最大的多次波检波波场能流密度矢量；
步骤七，通过步骤三获得的模值最大的多次波震源波场能流密度矢量与步骤六获得的模值最大的多次波检波波场能流密度矢量，计算地震波场的入射角；
步骤八，通过步骤一获得的多次波震源波场的位移场、步骤四获得的多次波检波波场的位移场以及所述步骤七获得的地震波场的入射角来计算多次波成像角度域共成像点道集。
作为本发明优选的方案，所述步骤一中，纵波介质速度场基于二阶声波波动方程，输入地表记录的炮集地震数据进行多次波地震波场的正向延拓，并保存每一空间点的不同传播时间的多次波震源波场位移场PS(x；t)，
采用如下二阶声波波动方程确定多次波位移场PS(x；t)：
式中，x＝x,z表示笛卡尔直角坐标系中一个空间点的坐标，x为水平方向的坐标，z为垂直方向的坐标；VP(x)为位于空间点x处的纵波速度；t为波场的传播时间；P为多次波波场的位移场，上标S表示震源波场；xR表示检波点的坐标，Dobs(xR；t)表示位于xR处的检波点在时间t时接收到的炮集地震数据。
作为本发明优选的方案，所述步骤二中，
采用如下多次波震源波场位移场梯度运算公式计算多次波震源波场位移场关于x方向的导数与多次波震源波场位移场关于z方向的导数：
式中，为在震源波场传播时间t位于空间点x处的多次波震源波场位移场关于x方向的导数，/>为在波场传播时间t位于空间点x处的多次波震源波场位移场关于z方向的导数；x＝x,z为笛卡尔直角坐标系中一个空间点的坐标，x为水平方向的坐标，z为垂直方向的坐标；t为震源波场的传播时间；P为多次波波场的位移场，上标S表示震源波场；
采用如下公式计算多次波震源波场位移场关于传播时间t的导数：
式中，为多次波震源波场位移场关于震源波场的传播时间t的导数，x为笛卡尔直角坐标系中一个空间点的坐标，t为震源波场的传播时间；P为多次波波场的位移场，上标S表示震源波场；
采用如下公式计算多次波震源波场的能流密度矢量：
式中，为多次波震源波场的能流密度矢量在x方向上的分量，/>为在震源波场传播时间t位于空间点x处的多次波震源波场位移场关于x方向的导数，/>为多次波震源波场位移场关于传播时间t的导数，x为笛卡尔直角坐标系中一个空间点的坐标；t为震源波场的传播时间；P为多次波波场的位移场，上标S表示震源波场；/>为多次波震源波场的能流密度矢量在z方向上的分量，/>为在震源波场传播时间t位于空间点x处的多次波震源波场位移场关于z方向的导数。
作为本发明优选的方案，所述步骤三中，
采用如下公式计算多次波扫描震源波场能流密度矢量：
式中，表示震源波场的传播时间，Tmax为最大传播时间，twin为给定的时间窗大小，/>表示震源波场传播/>时刻的多次波震源波场能流密度矢量在时间窗范围内的震源波场能流密度矢量模值最大的时刻；/>为计算的多次波震源波场能流密度矢量，/>表示多次波震源波场能流密度矢量在x方向上的分量，/>表示多次波震源波场能流密度矢量在z方向上的分量；/>为获得的多次波扫描震源波场能流密度矢量，/>表示多次波扫描震源波场能流密度矢量在x方向上的分量，/>表示多次波扫描震源波场能流密度矢量在z方向上的分量；上标S表示震源波场，||表示求模运算，max表示最大值，Δt为地震波场的传播时间的采样间隔；
根据上述公式首先计算震源波场传播时刻使得能流密度矢量最大时的时刻然后，将震源波场传播/>时刻的多次波扫描震源波场能流密度/>赋为计算的震源波场能流密度在传播时刻timax的值。
作为本发明优选的方案，所述步骤四中，纵波介质速度场基于二阶声波波动方程，输入地表记录的炮集地震数据进行多次波地震波场的反向延拓，并保存每一空间点的不同传播时间的多次波检波位移场PR(x；t)，
采用如下二阶声波波动方程确定多次波检波位移场PR(x；t)：
式中，x＝x,z表示笛卡尔直角坐标系中一个空间点的坐标，x为水平方向的坐标，z为垂直方向的坐标；VP(x)为位于空间点x处的纵波速度；t为波场的传播时间；P为多次波波场的位移场，上标R为检波波场；xR表示检波点的坐标，Dobs(xR；t)表示位于xR处的检波点在时间t时接收到的炮集地震数据。
作为本发明优选的方案，所述步骤五中，
采用如下多次波检波波场位移场梯度运算计算多次波检波波场位移场关于x方向的导数与多次波检波波场位移场关于z方向的导数，
式中，为在检波波场传播时间t位于空间点x处的多次波检波波场位移场关于x方向的导数；/>为在波场传播时间t位于空间点x处的多次波震源波场位移场关于z方向的导数；x＝x,z为笛卡尔直角坐标系中一个空间点的坐标，x为水平方向的坐标，z为垂直方向的坐标；t为检波波场的传播时间；P为多次波波场的位移场，上标R为检波波场；
采用如下公式计算多次波检波波场位移场关于传播时间t的导数：
式中，为多次波检波波场位移场关于检波波场的传播时间t的导数，x为笛卡尔直角坐标系中一个空间点的坐标，t为检波波场的传播时间，P为多次波波场的位移场，上标R为检波波场；
采用如下如下公式计算多次波检波波场的能流密度矢量：
式中，为多次波检波波场的能流密度矢量在x方向上的分量，/>为在检波波场传播时间t位于空间点x处的多次波检波波场位移场关于x方向的导数，/>为多次波检波波场位移场关于检波波场传播时间t的导数；x为笛卡尔直角坐标系中一个空间点的坐标；t为检波波场的传播时间；P为多次波波场的位移场，上标R为检波波场；/>为多次波检波波场的能流密度矢量在z方向上的分量，/>为在检波波场传播时间t位于空间点x处的多次波震源波场位移场关于z方向的导数。
作为本发明优选的方案，所述步骤六中，
采用如下公式计算多次波扫描检波波场能流密度矢量：
式中，表示检波波场的传播时间，Tmax为最大传播时间，twin为给定的时间窗大小，/>表示检波波场传播/>时刻的检波波场能流密度矢量在时间窗范围内的检波波场能流密度矢量模值最大的时刻，/>为计算的检波波场能流密度矢量，/>表示检波波场能流密度矢量在x方向上的分量，/>表示检波波场能流密度矢量在z方向上的分量；/>为获得的多次波扫描检波波场能流密度矢量，/>表示多次波扫描检波波场能流密度矢量在x方向上的分量，/>表示多次波扫描震源波场能流密度矢量在z方向上的分量；上标R表示检波波场，||表示求模运算，max表示求最大值运算；Δt为地震波场的传播时间的采样间隔；
根据上述公式首先计算检波波场传播时刻使得能流密度矢量最大时的时刻然后，将检波波场传播/>时刻的多次波扫描检波波场能流密度/>赋为计算的检波波场能流密度在传播时刻timax的值。
作为本发明优选的方案，所述步骤七中，
采用如下公式计算多次波波场在地下传播的入射角：：
式中，α(x；ti)表示地下空间点x在多次波波场传播时刻ti时的地震波场入射角；ti,i＝0～Tmax表示地震波场的传播时刻，Tmax为最大传播时间；表示震源波场的传播时间；/>表示检波波场的传播时间；/>为获得的多次波扫描震源波场能流密度矢量，/>表示多次波扫描震源波场能流密度矢量在x方向上的分量，表示多次波扫描震源波场能流密度矢量在z方向上的分量；/>为获得的多次波扫描检波波场能流密度矢量，/>表示多次波扫描检波波场能流密度矢量在x方向上的分量，/>表示多次波扫描震源波场能流密度矢量在z方向上的分量；上标S为震源波场，上标R为检波波场；·表示点积运算，||表示求模运算。
作为本发明优选的方案，所述步骤八中，
采用如下公式计算多次波成像的角度域共成像点道集：
式中，I表示多次波成像角度域共成像点道集，x为笛卡尔直角坐标系中一个空间点的坐标，αk为角度域共成像点道集中的入射角；t为波场的传播时间；P为多次波波场的位移场，上标S为震源波场，上标R为检波波场；*表示卷积运算符；e表示e指数函数；σ表示方差，取值为正有理数。
此外，本发明的另一方面还提供了多次波成像角度域共成像点道集的计算设备，其包括处理器、存储器，所述存储器用于存放至少一可执行指令，所述可执行指令使所述处理器执行所述的多次波成像角度域共成像点道集的计算方法对应的操作。
本发明实施例多次波成像角度域共成像点道集的计算方法及计算设备与现有技术相比，其有益效果在于：
本发明由于地震波场的能流密度矢量能够直接反映地震波场的传播方法，具有明确的物理意义，利用能流密度矢量计算多次波成像的角度域共成像点道集，保证多次波成像角度域共成像点道集的计算准确度高，能够提供地震勘探中多次波成像的角度域共成像点道集。
附图说明
为了更清楚地说明本发明实施例的技术方案，下面将对实施例的附图作简单地介绍。
图1是本发明实施例提供的步骤流程图；
图2a是本发明实施例提供的纵波介质速度场模型的示意图；
图2b是地表记录的炮集地震数据的示意图；
图3是正向延拓至1500ms时的多次波震源波场位移场的示意图；
图4a是正向延拓至1500ms时的多次波震源波场位移场关于x方向的导数的示意图；
图4b是正向延拓至1500ms时的多次波震源波场位移场关于z方向的导数的示意图；
图4c是正向延拓至1500ms时的多次波震源波场位移场关于传播时间的导数的示意图；
图4d是正向延拓至1500ms时的多次波震源波场的能流密度矢量在x方向上的分量的示意图；
图4e是正向延拓至1500ms时的多次波震源波场的能流密度矢量在z方向上的分量的示意图；
图5a是正向延拓至1500ms时的多次波扫描震源波场能流密度矢量在x方向上的分量的示意图；
图5b是正向延拓至1500ms时的多次波扫描震源波场能流密度矢量在z方向上的分量的示意图；
图6是反向延拓至1500ms时的多次波检波波场位移场的示意图；
图7a是反向延拓至1500ms时的多次波检波波场位移场关于x方向的导数的示意图；
图7b是反向延拓至1500ms时的多次波检波波场位移场关于z方向的导数的示意图；
图7c是反向延拓至1500ms时的多次波检波波场位移场关于传播时间的导数的示意图；
图7d是反向延拓至1500ms时的多次波检波波场能流密度矢量在x方向上的分量的示意图；
图7e是反向延拓至1500ms时的多次波检波波场能流密度矢量在z方向上的分量的示意图；
图8a是反向延拓至1500ms时的多次波扫描检波波场能流密度矢量在x方向上的分量的示意图；
图8b是反向延拓至1500ms时的多次波扫描检波波场能流密度矢量在z方向上的分量的示意图；
图9是反向延拓至1500ms时求得的各空间点的地震波场入射角的示意图；
图10是计算的水平位置位于图2a中1km处的多次波成像的角度域共成像点道集的示意图。
具体实施方式
下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。以下实施例用于说明本发明，但不用来限制本发明的范围。
在本发明的描述中，需要理解的是，术语“上”、“下”、“左”、“右”、“顶”、“底”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。应当理解的是，本发明中采用术语“第一”、“第二”等来描述各种信息，但这些信息不应限于这些术语，这些术语仅用来将同一类型的信息彼此区分开。例如，在不脱离本发明范围的情况下，“第一”信息也可以被称为“第二”信息，类似的，“第二”信息也可以被称为“第一”信息。
如图1至图10所示，本发明实施例优选实施例的多次波成像角度域共成像点道集的计算方法，包括以下步骤：包括以下步骤：
步骤一，基于声波波动方程，输入地表记录的炮集数据作为震源，进行地震波场的正向延拓获得多次波震源波场的位移场(S1)；在本实施例中，具体以图2a所示的纵波介质速度场模型为例，基于二阶声波波动方程，输入地表记录的炮集地震数据(如图2b所示)进行多次波地震波场的正向延拓，并保存每一空间点的不同传播时间的多次波震源波场位移场PS(x；t)，
采用如下二阶声波波动方程确定多次波位移场PS(x；t)：
式中，x＝x,z表示笛卡尔直角坐标系中一个空间点的坐标，x为水平方向的坐标，z为垂直方向的坐标；VP(x)为位于空间点x处的纵波速度；t为波场的传播时间；P为多次波波场的位移场，上标S表示震源波场；xR表示检波点的坐标，Dobs(xR；t)表示位于xR处的检波点在时间t时接收到的炮集地震数据。这样能够获得正向延拓至1500ms时如图3所示的多次波震源波场位移场。
步骤二，将步骤一保存的多次波震源波场位移场进行梯度运算和关于传播时间的导数运算，并根据能流密度的计算公式获得震源波场的能流密度矢量(S2)；具体的，在所述步骤二中，
采用如下多次波震源波场位移场梯度运算公式计算多次波震源波场位移场关于x方向的导数与多次波震源波场位移场关于z方向的导数：
式中，为在震源波场传播时间t位于空间点x处的多次波震源波场位移场关于x方向的导数，/>为在波场传播时间t位于空间点x处的多次波震源波场位移场关于z方向的导数；x＝x,z为笛卡尔直角坐标系中一个空间点的坐标，x为水平方向的坐标，z为垂直方向的坐标；t为震源波场的传播时间；P为多次波波场的位移场，上标S表示震源波场；
在本实施例中，这样能够获得正向延拓至1500ms时如图4a所示的多次波震源波场位移场关于x方向的导数，以及，正向延拓至1500ms时如图4b所示的多次波震源波场位移场关于z方向的导数；
采用如下公式计算多次波震源波场位移场关于传播时间t的导数：
式中，为多次波震源波场位移场关于震源波场的传播时间t的导数，x为笛卡尔直角坐标系中一个空间点的坐标，t为震源波场的传播时间；P为多次波波场的位移场，上标S表示震源波场；
在本实施例中，这样能够获得正向延拓至1500ms时如图4c所示的多次波震源波场位移场关于传播时间t的导数；
采用如下公式计算多次波震源波场的能流密度矢量：
式中，为多次波震源波场的能流密度矢量在x方向上的分量，/>为在震源波场传播时间t位于空间点x处的多次波震源波场位移场关于x方向的导数，/>为多次波震源波场位移场关于传播时间t的导数，x为笛卡尔直角坐标系中一个空间点的坐标；t为震源波场的传播时间；P为多次波波场的位移场，上标S表示震源波场；/>为多次波震源波场的能流密度矢量在z方向上的分量，/>为在震源波场传播时间t位于空间点x处的多次波震源波场位移场关于z方向的导数；
在本实施例中，这样能够获得正向延拓至1500ms时如图4d所示的多次波震源波场的能流密度矢量在x方向上的分量，以及，获得正向延拓至1500ms时如图4e所示的多次波震源波场的能流密度矢量在z方向上的分量。
步骤三，在给定的时间窗内，扫描每一时刻步骤二获得的震源波场能流密度矢量的模值，并保存模值最大的震源波场能流密度矢量(S3)；具体的，在所述步骤三中，
采用如下公式计算多次波扫描震源波场能流密度矢量：
式中，表示震源波场的传播时间，Tmax为最大传播时间，twin为给定的时间窗大小，/>表示震源波场传播/>时刻的多次波震源波场能流密度矢量在时间窗范围内的震源波场能流密度矢量模值最大的时刻；/>为计算的多次波震源波场能流密度矢量，/>表示多次波震源波场能流密度矢量在x方向上的分量，/>表示多次波震源波场能流密度矢量在z方向上的分量；/>为获得的多次波扫描震源波场能流密度矢量，/>表示多次波扫描震源波场能流密度矢量在x方向上的分量，/>表示多次波扫描震源波场能流密度矢量在z方向上的分量；上标S表示震源波场，||表示求模运算，max表示最大值，Δt为地震波场的传播时间的采样间隔；
根据上述公式首先计算震源波场传播时刻使得能流密度矢量最大时的时刻然后，将震源波场传播/>时刻的多次波扫描震源波场能流密度/>赋为计算的震源波场能流密度在传播时刻timax的值。在本实施例中，这样能够获得正向延拓至1500ms时如图5a所示的多次波扫描震源波场能流密度矢量在x方向上的分量，获得正向延拓至1500ms时如图5b所示的多次波扫描震源波场能流密度矢量在z方向上的分量。
步骤四，基于声波波动方程，输入地表记录的炮集数据作为边值条件，进行地震波场的逆向延拓获得多次波检波波场的位移场(S4)；具体的，在所述步骤四中，纵波介质速度场(如图2a所示)基于二阶声波波动方程，输入地表记录的炮集地震数据(如图2b所示)进行多次波地震波场的反向延拓，并保存每一空间点的不同传播时间的多次波检波位移场PR(x；t)，
采用如下二阶声波波动方程确定多次波检波位移场PR(x；t)：
式中，x＝x,z表示笛卡尔直角坐标系中一个空间点的坐标，x为水平方向的坐标，z为垂直方向的坐标；VP(x)为位于空间点x处的纵波速度；t为波场的传播时间；P为多次波波场的位移场，上标R为检波波场；xR表示检波点的坐标，Dobs(xR；t)表示位于xR处的检波点在时间t时接收到的炮集地震数据。在本实施例中，这样能够获得反向延拓至1500ms时如图6所示的多次波检波波场位移场。
步骤五，将步骤四保存的多次波检波波场位移场进行梯度运算和关于传播时间的导数运算，并根据能流密度的计算公式获得多次波检波波场的能流密度矢量(S5)；具体的，在所述步骤五中，采用如下多次波检波波场位移场梯度运算计算多次波检波波场位移场关于x方向的导数与多次波检波波场位移场关于z方向的导数，
式中，为在检波波场传播时间t位于空间点x处的多次波检波波场位移场关于x方向的导数；/>为在波场传播时间t位于空间点x处的多次波震源波场位移场关于z方向的导数；x＝x,z为笛卡尔直角坐标系中一个空间点的坐标，x为水平方向的坐标，z为垂直方向的坐标；t为检波波场的传播时间；P为多次波波场的位移场，上标R为检波波场；
在本实施例中，这样能够获得反向延拓至1500ms时如图7a所示的多次波检波波场位移场关于x方向的导数，以及反向延拓至1500ms时如图7b所示的多次波检波波场位移场关于z方向的导数；
采用如下公式计算多次波检波波场位移场关于传播时间t的导数：
式中，为多次波检波波场位移场关于检波波场的传播时间t的导数，x为笛卡尔直角坐标系中一个空间点的坐标，t为检波波场的传播时间，P为多次波波场的位移场，上标R为检波波场；
在本实施例中，这样能够获得反向延拓至1500ms时如图7c所示的多次波检波波场位移场关于传播时间t的导数；
采用如下如下公式计算多次波检波波场的能流密度矢量：
式中，为多次波检波波场的能流密度矢量在x方向上的分量，/>为在检波波场传播时间t位于空间点x处的多次波检波波场位移场关于x方向的导数，/>为多次波检波波场位移场关于检波波场传播时间t的导数；x为笛卡尔直角坐标系中一个空间点的坐标；t为检波波场的传播时间；P为多次波波场的位移场，上标R为检波波场；/>为多次波检波波场的能流密度矢量在z方向上的分量，/>为在检波波场传播时间t位于空间点x处的多次波震源波场位移场关于z方向的导数；
在本实施例中，这样能够获得反向延拓至1500ms时如图7d所示的多次波检波波场能流密度矢量在x方向上的分量，以及，反向延拓至1500ms时如图7e所示的多次波检波波场能流密度矢量在z方向上的分量。
步骤六，在给定的时间窗内，扫描每一时刻步骤五获得的检波波场能流密度矢量的模值，并保存模值最大的多次波检波波场能流密度矢量(S6)；具体的，在所述步骤六中，
采用如下公式计算多次波扫描检波波场能流密度矢量：
式中，表示检波波场的传播时间，Tmax为最大传播时间，twin为给定的时间窗大小，/>表示检波波场传播/>时刻的检波波场能流密度矢量在时间窗范围内的检波波场能流密度矢量模值最大的时刻，/>为计算的检波波场能流密度矢量，/>表示检波波场能流密度矢量在x方向上的分量，/>表示检波波场能流密度矢量在z方向上的分量；/>为获得的多次波扫描检波波场能流密度矢量，/>表示多次波扫描检波波场能流密度矢量在x方向上的分量，/>表示多次波扫描震源波场能流密度矢量在z方向上的分量；上标R表示检波波场，||表示求模运算，max表示求最大值运算；Δt为地震波场的传播时间的采样间隔；
根据上述公式首先计算检波波场传播时刻使得能流密度矢量最大时的时刻然后，将检波波场传播/>时刻的多次波扫描检波波场能流密度/>赋为计算的检波波场能流密度在传播时刻timax的值。在本实施例中，这样能够获得反向延拓至1500ms时如图8a所示的多次波扫描检波波场能流密度矢量在x方向上的分量，以及，反向延拓至1500ms时如图8b所示的多次波扫描检波波场能流密度矢量在z方向上的分量。
步骤七，通过步骤三获得的模值最大的多次波震源波场能流密度矢量与步骤六获得的模值最大的多次波检波波场能流密度矢量，计算地震波场的入射角(S7)；具体的，在所述步骤七中，
采用如下公式计算多次波波场在地下传播的入射角：
式中，α(x；ti)表示地下空间点x在多次波波场传播时刻ti时的地震波场入射角；ti,i＝0～Tmax表示地震波场的传播时刻，Tmax为最大传播时间；表示震源波场的传播时间；/>表示检波波场的传播时间；/>为获得的多次波扫描震源波场能流密度矢量，/>表示多次波扫描震源波场能流密度矢量在x方向上的分量，表示多次波扫描震源波场能流密度矢量在z方向上的分量；/>为获得的多次波扫描检波波场能流密度矢量，/>表示多次波扫描检波波场能流密度矢量在x方向上的分量，/>表示多次波扫描震源波场能流密度矢量在z方向上的分量；上标S为震源波场，上标R为检波波场；·表示点积运算，||表示求模运算；在本实施例中，这样能够获得反向延拓至1500ms时求得的各空间点的地震波场入射角(如图9所示)。
步骤八，通过步骤一获得的多次波震源波场的位移场、步骤四获得的多次波检波波场的位移场以及所述步骤七获得的地震波场的入射角来计算多次波成像角度域共成像点道集(S8)；具体的，所述步骤八中，
采用如下公式计算多次波成像的角度域共成像点道集：
式中，I表示多次波成像角度域共成像点道集，x为笛卡尔直角坐标系中一个空间点的坐标，αk为角度域共成像点道集中的入射角；t为波场的传播时间；P为多次波波场的位移场，上标S为震源波场，上标R为检波波场；*表示卷积运算符；e表示e指数函数；σ表示方差，取值为正有理数。在本实施例中，如图10所示为计算的水平位置位于图2a中1km处的多次波成像的角度域共成像点道集。
综上，本发明由于地震波场的能流密度矢量能够直接反映地震波场的传播方法，具有明确的物理意义，利用能流密度矢量计算多次波成像的角度域共成像点道集，保证多次波成像角度域共成像点道集的计算准确度高，能够提供地震勘探中多次波成像的角度域共成像点道集。
此外，本发明还提供了地震勘探中的多次波成像角度域共成像点道集的计算设备，其包括处理器、存储器，所述存储器用于存放至少一可执行指令，所述可执行指令使所述处理器执行所述的多次波成像角度域共成像点道集的计算方法对应的操作，具体的操作步骤为，
步骤一，基于声波波动方程，输入地表记录的炮集数据作为震源，进行地震波场的正向延拓获得多次波震源波场的位移场；
步骤二，将步骤一保存的多次波震源波场位移场进行梯度运算和关于传播时间的导数运算，并根据能流密度的计算公式获得震源波场的能流密度矢量；
步骤三，在给定的时间窗内，扫描每一时刻步骤二获得的震源波场能流密度矢量的模值，并保存模值最大的震源波场能流密度矢量；
步骤四，基于声波波动方程，输入地表记录的炮集数据作为边值条件，进行地震波场的逆向延拓获得多次波检波波场的位移场；
步骤五，将步骤四保存的多次波检波波场位移场进行梯度运算和关于传播时间的导数运算，并根据能流密度的计算公式获得多次波检波波场的能流密度矢量；
步骤六，在给定的时间窗内，扫描每一时刻步骤五获得的检波波场能流密度矢量的模值，并保存模值最大的多次波检波波场能流密度矢量；
步骤七，通过步骤三获得的模值最大的多次波震源波场能流密度矢量与步骤六获得的模值最大的多次波检波波场能流密度矢量，计算地震波场的入射角；
步骤八，通过步骤一获得的多次波震源波场的位移场、步骤四获得的多次波检波波场的位移场以及所述步骤七获得的地震波场的入射角来计算多次波成像角度域共成像点道集。
另外，本发明还提供了地震勘探中的多次波成像角度域共成像点道集的计算装置，包括：
多次波震源波场位移场获取模块，基于声波波动方程，输入地表记录的炮集数据作为震源，进行地震波场的正向延拓获得多次波震源波场的位移场；
震源波场能流密度矢量获取模块，将保存的多次波震源波场位移场进行梯度运算和关于传播时间的导数运算，并根据能流密度的计算公式获得震源波场的能流密度矢量；
扫描震源波场能流密度矢量获取模块，在给定的时间窗内，扫描每一时刻震源波场能流密度矢量的模值，并保存模值最大的震源波场能流密度矢量；
多次波检波波场位移场获取模块，基于声波波动方程，输入地表记录的炮集数据作为边值条件，进行地震波场的逆向延拓获得多次波检波波场的位移场；
检波波场能流密度矢量获取模块，将保存的多次波检波波场位移场进行梯度运算和关于传播时间的导数运算，并根据能流密度的计算公式获得多次波检波波场的能流密度矢量；
扫描检波波场能流密度矢量获取模块，在给定的时间窗内，扫描每一时刻步骤五获得的检波波场能流密度矢量的模值，并保存模值最大的多次波检波波场能流密度矢量；
地震波场入射角获取模块，通过模值最大的多次波震源波场能流密度矢量与模值最大的多次波检波波场能流密度矢量，计算地震波场的入射角；
多次波角度域共成像点道集获取模块；通过多次波震源波场的位移场、多次波检波波场的位移场以及地震波场的入射角来计算多次波成像角度域共成像点道集。
在本发明的描述中，需要说明的是，除非另有明确的规定和限定，术语“相连”、“连接”应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通。对于本领域的普通技术人员而言，可以具体情况理解上述术语在本发明中的具体含义。
以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和替换，这些改进和替换也应视为本发明的保护范围。

Claims (10)

1.多次波成像角度域共成像点道集的计算方法，其特征在于，包括以下步骤：
步骤一，基于二阶声波波动方程，输入地表记录的炮集数据作为震源，进行地震波场的正向延拓获得多次波震源波场的位移场；
步骤二，将步骤一保存的多次波震源波场位移场进行梯度运算和关于传播时间的导数运算，并根据能流密度的计算公式获得震源波场的能流密度矢量；
步骤三，在给定的时间窗内，扫描每一时刻步骤二获得的震源波场能流密度矢量的模值，并保存模值最大的震源波场能流密度矢量；
步骤四，基于二阶声波波动方程，输入地表记录的炮集数据作为边值条件，进行地震波场的逆向延拓获得多次波检波波场的位移场；
步骤五，将步骤四保存的多次波检波波场位移场进行梯度运算和关于传播时间的导数运算，并根据能流密度的计算公式获得多次波检波波场的能流密度矢量；
步骤六，在给定的时间窗内，扫描每一时刻步骤五获得的检波波场能流密度矢量的模值，并保存模值最大的多次波检波波场能流密度矢量；
步骤七，通过步骤三获得的模值最大的多次波震源波场能流密度矢量与步骤六获得的模值最大的多次波检波波场能流密度矢量，计算地震波场的入射角；
步骤八，通过步骤一获得的多次波震源波场的位移场、步骤四获得的多次波检波波场的位移场以及所述步骤七获得的地震波场的入射角来计算多次波成像角度域共成像点道集。
2.如权利要求1所述的多次波成像角度域共成像点道集的计算方法，其特征在于，所述步骤一中，纵波介质速度场基于二阶声波波动方程，输入地表记录的炮集地震数据进行多次波地震波场的正向延拓，并保存每一空间点的不同传播时间的多次波震源波场位移场PS(x；t)，
采用如下二阶声波波动方程确定多次波位移场PS(x；t)：
式中，x＝x,z表示笛卡尔直角坐标系中一个空间点的坐标，x为水平方向的坐标，z为垂直方向的坐标；VP(x)为位于空间点x处的纵波速度；t为波场的传播时间；P为多次波波场的位移场，上标S表示震源波场；xR表示检波点的坐标，Dobs(xR；t)表示位于xR处的检波点在时间t时接收到的炮集地震数据。
3.如权利要求2所述的多次波成像角度域共成像点道集的计算方法，其特征在于，所述步骤二中，
采用如下多次波震源波场位移场梯度运算公式计算多次波震源波场位移场关于x方向的导数与多次波震源波场位移场关于z方向的导数：
式中，为在震源波场传播时间t位于空间点x处的多次波震源波场位移场关于x方向的导数，/>为在波场传播时间t位于空间点x处的多次波震源波场位移场关于z方向的导数；x＝x,z为笛卡尔直角坐标系中一个空间点的坐标，x为水平方向的坐标，z为垂直方向的坐标；t为震源波场的传播时间；P为多次波波场的位移场，上标S表示震源波场；
采用如下公式计算多次波震源波场位移场关于传播时间t的导数：
式中，为多次波震源波场位移场关于震源波场的传播时间t的导数，x为笛卡尔直角坐标系中一个空间点的坐标，t为震源波场的传播时间；P为多次波波场的位移场，上标S表示震源波场；
采用如下公式计算多次波震源波场的能流密度矢量：
式中，为多次波震源波场的能流密度矢量在x方向上的分量，/>为在震源波场传播时间t位于空间点x处的多次波震源波场位移场关于x方向的导数，/>为多次波震源波场位移场关于传播时间t的导数，x为笛卡尔直角坐标系中一个空间点的坐标；t为震源波场的传播时间；P为多次波波场的位移场，上标S表示震源波场；/>为多次波震源波场的能流密度矢量在z方向上的分量，/>为在震源波场传播时间t位于空间点x处的多次波震源波场位移场关于z方向的导数。
4.如权利要求3所述的多次波成像角度域共成像点道集的计算方法，其特征在于，所述步骤三中，
采用如下公式计算多次波扫描震源波场能流密度矢量：
式中，表示震源波场的传播时间，Tmax为最大传播时间，twin为给定的时间窗大小，/>表示震源波场传播/>时刻的多次波震源波场能流密度矢量在时间窗范围内的震源波场能流密度矢量模值最大的时刻；/>为计算的多次波震源波场能流密度矢量，/>表示多次波震源波场能流密度矢量在x方向上的分量，/>表示多次波震源波场能流密度矢量在z方向上的分量；/>为获得的多次波扫描震源波场能流密度矢量，/>表示多次波扫描震源波场能流密度矢量在x方向上的分量，/>表示多次波扫描震源波场能流密度矢量在z方向上的分量；上标S表示震源波场，| |表示求模运算，max表示最大值，Δt为地震波场的传播时间的采样间隔；
根据上述公式首先计算震源波场传播时刻使得能流密度矢量最大时的时刻/>然后，将震源波场传播/>时刻的多次波扫描震源波场能流密度/>赋为计算的震源波场能流密度在传播时刻timax的值。
5.如权利要求4所述的多次波成像角度域共成像点道集的计算方法，其特征在于，所述步骤四中，纵波介质速度场基于二阶声波波动方程，输入地表记录的炮集地震数据进行多次波地震波场的反向延拓，并保存每一空间点的不同传播时间的多次波检波位移场PR(x；t)，
采用如下二阶声波波动方程确定多次波检波位移场PR(x；t)：
式中，x＝x,z表示笛卡尔直角坐标系中一个空间点的坐标，x为水平方向的坐标，z为垂直方向的坐标；VP(x)为位于空间点x处的纵波速度；t为波场的传播时间；P为多次波波场的位移场，上标R为检波波场；xR表示检波点的坐标，Dobs(xR；t)表示位于xR处的检波点在时间t时接收到的炮集地震数据。
6.如权利要求5所述的多次波成像角度域共成像点道集的计算方法，其特征在于，所述步骤五中，
采用如下多次波检波波场位移场梯度运算计算多次波检波波场位移场关于x方向的导数与多次波检波波场位移场关于z方向的导数，
式中，为在检波波场传播时间t位于空间点x处的多次波检波波场位移场关于x方向的导数；/>为在波场传播时间t位于空间点x处的多次波震源波场位移场关于z方向的导数；x为笛卡尔直角坐标系中一个空间点的坐标，x为水平方向的坐标，z为垂直方向的坐标；t为检波波场的传播时间；P为多次波波场的位移场，上标R为检波波场；
采用如下公式计算多次波检波波场位移场关于传播时间t的导数：
式中，为多次波检波波场位移场关于检波波场的传播时间t的导数，x为笛卡尔直角坐标系中一个空间点的坐标，t为检波波场的传播时间，P为多次波波场的位移场，上标R为检波波场；
采用如下公式计算多次波检波波场的能流密度矢量：
式中，为多次波检波波场的能流密度矢量在x方向上的分量，/>为在检波波场传播时间t位于空间点x处的多次波检波波场位移场关于x方向的导数，/>为多次波检波波场位移场关于检波波场传播时间t的导数；x为笛卡尔直角坐标系中一个空间点的坐标；t为检波波场的传播时间；P为多次波波场的位移场，上标R为检波波场；/>为多次波检波波场的能流密度矢量在z方向上的分量，/>为在检波波场传播时间t位于空间点x处的多次波震源波场位移场关于z方向的导数。
7.如权利要求6所述的多次波成像角度域共成像点道集的计算方法，其特征在于，所述步骤六中，
采用如下公式计算多次波扫描检波波场能流密度矢量：
式中，表示检波波场的传播时间，Tmax为最大传播时间，twin为给定的时间窗大小，/>表示检波波场传播/>时刻的检波波场能流密度矢量在时间窗范围内的检波波场能流密度矢量模值最大的时刻，/>为计算的检波波场能流密度矢量，/>表示检波波场能流密度矢量在x方向上的分量，/>表示检波波场能流密度矢量在z方向上的分量；/>为获得的多次波扫描检波波场能流密度矢量，/>表示多次波扫描检波波场能流密度矢量在x方向上的分量，/>表示多次波扫描震源波场能流密度矢量在z方向上的分量；上标R表示检波波场，| |表示求模运算，max表示求最大值运算；Δt为地震波场的传播时间的采样间隔；
根据上述公式首先计算检波波场传播时刻使得能流密度矢量最大时的时刻/>然后，将检波波场传播/>时刻的多次波扫描检波波场能流密度/>赋为计算的检波波场能流密度在传播时刻timax的值。
8.如权利要求7所述的多次波成像角度域共成像点道集的计算方法，其特征在于，所述步骤七中，
采用如下公式计算多次波波场在地下传播的入射角：
式中，α(x；ti)表示地下空间点x在多次波波场传播时刻ti时的地震波场入射角；ti,i＝0～Tmax表示地震波场的传播时刻，Tmax为最大传播时间；表示震源波场的传播时间；/>表示检波波场的传播时间；/>为获得的多次波扫描震源波场能流密度矢量，/>表示多次波扫描震源波场能流密度矢量在x方向上的分量，/>表示多次波扫描震源波场能流密度矢量在z方向上的分量；/>为获得的多次波扫描检波波场能流密度矢量，/>表示多次波扫描检波波场能流密度矢量在x方向上的分量，表示多次波扫描震源波场能流密度矢量在z方向上的分量；上标S为震源波场，上标R为检波波场；·表示点积运算，| |表示求模运算。
9.如权利要求8所述的多次波成像角度域共成像点道集的计算方法，其特征在于，所述步骤八中，
采用如下公式计算多次波成像的角度域共成像点道集：
式中，I表示多次波成像角度域共成像点道集，x为笛卡尔直角坐标系中一个空间点的坐标，αk为角度域共成像点道集中的入射角；t为波场的传播时间；P为多次波波场的位移场，上标S为震源波场，上标R为检波波场；*表示卷积运算符；e表示e指数函数；σ表示方差，取值为正有理数。
10.地震勘探中的多次波成像角度域共成像点道集的计算设备，其特征在于，包括处理器、存储器，所述存储器用于存放至少一可执行指令，所述可执行指令使所述处理器执行如权利要求1～9中任一项所述的多次波成像角度域共成像点道集的计算方法对应的操作。
CN202211121320.3A 2022-09-15 2022-09-15 多次波成像角度域共成像点道集计算方法及计算设备 Active CN115951407B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211121320.3A CN115951407B (zh) 2022-09-15 2022-09-15 多次波成像角度域共成像点道集计算方法及计算设备

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211121320.3A CN115951407B (zh) 2022-09-15 2022-09-15 多次波成像角度域共成像点道集计算方法及计算设备

Publications (2)

Publication Number Publication Date
CN115951407A CN115951407A (zh) 2023-04-11
CN115951407B true CN115951407B (zh) 2023-09-29

Family

ID=87288316

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211121320.3A Active CN115951407B (zh) 2022-09-15 2022-09-15 多次波成像角度域共成像点道集计算方法及计算设备

Country Status (1)

Country Link
CN (1) CN115951407B (zh)

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102841375A (zh) * 2012-09-06 2012-12-26 中国石油大学(华东) 一种复杂条件下基于角度域共成像点道集的层析速度反演方法
CN103616721A (zh) * 2013-11-25 2014-03-05 中国石油天然气股份有限公司 基于二阶偏微分波动方程的pml吸收边界条件的方法
CN107894613A (zh) * 2017-10-26 2018-04-10 中国石油天然气集团公司 弹性波矢量成像方法、装置、存储介质及设备
CN114428341A (zh) * 2020-09-30 2022-05-03 中国石油化工股份有限公司 振幅保真地震成像方法、装置、电子设备及介质
CN114647003A (zh) * 2022-03-18 2022-06-21 中山大学 转换波角度域共成像点道集剩余时差计算方法和设备

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2011152928A1 (en) * 2010-06-02 2011-12-08 Exxonmobil Upstream Research Company Efficient computation of wave equation migration angle gathers
US10557954B2 (en) * 2017-06-12 2020-02-11 Saudi Arabian Oil Company Modeling angle domain common image gathers from reverse time migration

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102841375A (zh) * 2012-09-06 2012-12-26 中国石油大学(华东) 一种复杂条件下基于角度域共成像点道集的层析速度反演方法
CN103616721A (zh) * 2013-11-25 2014-03-05 中国石油天然气股份有限公司 基于二阶偏微分波动方程的pml吸收边界条件的方法
CN107894613A (zh) * 2017-10-26 2018-04-10 中国石油天然气集团公司 弹性波矢量成像方法、装置、存储介质及设备
CN114428341A (zh) * 2020-09-30 2022-05-03 中国石油化工股份有限公司 振幅保真地震成像方法、装置、电子设备及介质
CN114647003A (zh) * 2022-03-18 2022-06-21 中山大学 转换波角度域共成像点道集剩余时差计算方法和设备

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
波动方程偏移角度域共成像道集计算方法;Sava P.;唐祥功;;油气地球物理(第02期);75-79 *

Also Published As

Publication number Publication date
CN115951407A (zh) 2023-04-11

Similar Documents

Publication Publication Date Title
US11092707B2 (en) Determining a component of a wave field
US8396668B2 (en) Marine seismic surveying employing interpolated multicomponent streamer pressure data
Schuster et al. A theoretical overview of model-based and correlation-based redatuming methods
Slob et al. Seismic reflector imaging using internal multiples with Marchenko-type equations
US8559266B2 (en) Seismic processing for the elimination of multiple reflections
AU2007272702B2 (en) Time lapse marine seismic surveying employing interpolated multicomponent streamer pressure data
RU2616650C2 (ru) Способ и устройство для обработки сейсмических данных
EP1879052A2 (en) Time lapse marine seismic surveying employing interpolated multicomponent streamer pressure data
US7460437B2 (en) Seismic data processing method and system for migration of seismic signals incorporating azimuthal variations in the velocity
EA032186B1 (ru) Сейсмическая адаптивная фокусировка
EA022531B1 (ru) Способ определения сейсмического атрибута по сейсмическим сигналам
SG195507A1 (en) Surface-related multiple elimination for depth-varying streamer
Li et al. High-resolution adaptive beamforming for borehole acoustic reflection imaging
CN115951407B (zh) 多次波成像角度域共成像点道集计算方法及计算设备
CN115469362B (zh) 一种地震勘探中的能流密度矢量计算方法
US7069200B2 (en) Method intended to obtain reflection travel times from an interpretation of migrated cylindrical wave seismic data
US6418380B1 (en) Method for seismic processing and in particular for three-dimensional seismic exploration using seismic data migration
WO2019118850A1 (en) Subsalt imaging tool for interpreters
Tsvankin et al. 3D description and inversion of reflection moveout of PS‐waves in anisotropic media
Yang et al. Least-squares reverse time migration with velocity errors
CN117368969A (zh) 一种基于近场记录模拟远场子波的方法、装置和设备
WO2019145746A1 (en) Method for generating a reconstructed seismic signal
EP3092520A1 (en) Determining a component of a wave field

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