CN117075197B - 一种横向各向同性介质波场分离的横波解耦方程构建方法 - Google Patents
一种横向各向同性介质波场分离的横波解耦方程构建方法 Download PDFInfo
- Publication number
- CN117075197B CN117075197B CN202311316458.3A CN202311316458A CN117075197B CN 117075197 B CN117075197 B CN 117075197B CN 202311316458 A CN202311316458 A CN 202311316458A CN 117075197 B CN117075197 B CN 117075197B
- Authority
- CN
- China
- Prior art keywords
- wave
- transverse
- stress
- decoupling
- velocity
- 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
- 238000000926 separation method Methods 0.000 title claims abstract description 40
- 238000010276 construction Methods 0.000 title claims abstract description 13
- 238000000034 method Methods 0.000 claims abstract description 96
- 239000002245 particle Substances 0.000 claims abstract description 70
- 239000011159 matrix material Substances 0.000 claims abstract description 55
- 230000008569 process Effects 0.000 claims description 28
- 238000005070 sampling Methods 0.000 claims description 16
- 229910000831 Steel Inorganic materials 0.000 claims description 10
- 239000010959 steel Substances 0.000 claims description 10
- 238000004364 calculation method Methods 0.000 claims description 5
- 238000000354 decomposition reaction Methods 0.000 claims description 3
- 239000006185 dispersion Substances 0.000 claims description 3
- 238000003384 imaging method Methods 0.000 description 10
- 239000010410 layer Substances 0.000 description 8
- 230000008878 coupling Effects 0.000 description 4
- 238000010168 coupling process Methods 0.000 description 4
- 238000005859 coupling reaction Methods 0.000 description 4
- 238000010586 diagram Methods 0.000 description 4
- 238000004458 analytical method Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000004880 explosion Methods 0.000 description 2
- 230000010287 polarization Effects 0.000 description 2
- 230000006978 adaptation Effects 0.000 description 1
- 238000007792 addition Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000004422 calculation algorithm Methods 0.000 description 1
- 238000012512 characterization method Methods 0.000 description 1
- 239000002355 dual-layer Substances 0.000 description 1
- 239000010910 field residue Substances 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000035939 shock Effects 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
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
-
- 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
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
-
- 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
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Remote Sensing (AREA)
- Life Sciences & Earth Sciences (AREA)
- Mathematical Physics (AREA)
- Computational Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Geology (AREA)
- Geophysics (AREA)
- Environmental & Geological Engineering (AREA)
- General Life Sciences & Earth Sciences (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Acoustics & Sound (AREA)
- Pure & Applied Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computing Systems (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Algebra (AREA)
- Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
Abstract
本发明提供了一种横向各向同性介质波场分离的横波解耦方程构建方法，涉及地球物理学技术领域。本发明基于弹性波一阶速度‑应力方程建立交错网格高阶有限差分格式，先对横向各向同性介质进行弹性波场延拓，再基于方程解耦将弹性刚度系数近似分解为待定的纵横波刚度矩阵弹性参数，引入修正伪亥姆霍兹算子，基于横波无散的反证法求解待定的横波刚度矩阵弹性参数，利用横波刚度矩阵弹性参数构造qS波应力并求取qS波质点振动速度，构造出qS波解耦一阶速度‑应力方程后，将波场延拓的质点振动速度场代入qS波和qP波的解耦一阶速度‑应力方程中，解耦分离纵横波波场，在时空域上实现了横向各向同性介质中弹性波场的无串扰分离。
Description
技术领域
本发明涉及地球物理学技术领域，具体涉及一种横向各向同性介质波场分离的横波解耦方程构建方法。
背景技术
地下介质普遍具有各向异性，考虑地球介质的各向异性特征有助于降低井震深度误差，为后续储层精细表征提供可靠的成像数据。作为最典型的各向异性介质类型的横向各向同性（VTI）介质，影响其弹性波成像精度的纵横波波场分离方法是目前的研究热点之一。
相较于各向同性介质，各向异性地震波场的复杂性与模型参数的增加为弹性波场分离的研究带来许多挑战，现有针对各向异性介质的波场分离理论尚不完善，分离算法复杂且计算效率较低。目前，各向异性介质下的波场分离方法主要分为两大类：第一类方法是波数域分离方法；这类方法根据纵横波的偏振方向，将波场投影到其偏振方向上从而实现纵横波解耦，但这类方法通常需要大量的傅里叶变换，计算量非常大，并且波数域分离方法通常应用于理论分析，无法用于常规时空域的波场模拟，在生产应用中鲜有研究。第二类方法是时空域分离方法；该类方法基于横向各向同性介质中解耦的纵横波方程，在时空域波场延拓的过程中实现波场分离，通常是给出纵波解耦方程，利用横波波场通过总波场减去纵波波场求差得到，这类方法是在时空域中进行的，得到了广泛的应用，但是，由于横波波场是通过求差得到的，横波波场之中呈现出纵横波耦合波场的残留，造成分离得到的横波波场存在严重串扰。
因此，针对横向各向同性介质弹性波波场分离中横波波场耦合严重的问题，亟需提出一种横向各向同性介质波场分离的横波解耦方程构建方法，实现时空域中对横向各向同性介质弹性波场的无串扰分离，为多分量地震数据的各向异性弹性波成像与反演奠定基础。
发明内容
本发明旨在解决横向各向同性介质弹性波波场分离中横波波场耦合严重的问题，提出了一种横向各向同性介质波场分离的横波解耦方程构建方法，基于横波无散的反证法在时空域上构建横向各向同性介质的qS波解耦一阶速度-应力方程，结合qP波解耦一阶速度-应力方程，实现了时空域上对横向各向同性介质中弹性波场的无串扰分离，为多分量地震数据的各向异性弹性波成像与反演奠定了基础。
为了实现上述目的，本发明采用如下技术方案：
一种横向各向同性介质波场分离的横波解耦方程构建方法，具体包括以下步骤：
步骤1，输入横向各向同性介质模型，对横向各向同性介质进行弹性波场延拓；
步骤2，基于方程解耦将弹性刚度系数近似分解为待定的纵横波刚度矩阵弹性参数；
步骤3，引入修正伪亥姆霍兹算子，基于横波无散的反证法求解待定的横波刚度矩阵弹性参数，利用横波刚度矩阵弹性参数构造qS波应力，根据qS波应力求取qS波质点振动速度，构造qS波解耦一阶速度-应力方程；
步骤4，将波场延拓得到的质点振动速度场代入qS波解耦一阶速度-应力方程和qP波解耦一阶速度-应力方程中，解耦分离纵横波波场。
优选地，所述步骤1中，将横向各向同性介质弹性波一阶速度-应力方程表示为：
（1）
（2）
式中，为横向各向同性介质模型的横向方向，为横向各向同性介质模型的纵
向方向，为横向各向同性介质模型的垂直方向；为横向各向同性介质的密度；
为指标索引，取值为；代表质点振动速度，取值为,其中，为质
点振动速度矢量场沿方向的分量，为质点振动速度矢量场沿方向的分量，为质点
振动速度矢量场沿方向的分量；为指标索引所对应质点振动速度在时间上的一阶导
数；为指标索引所对应的质点振动速度对方向求一阶导数，为指标索引所对
应的质点振动速度对方向求一阶导数，代表应力，取值为，其中，、、均为正应力，、、均为切
应力；为指标索引所对应的应力对方向求一阶导数；为指标索引所对应的应
力对时间求一阶导数；为横向各向同性介质刚度矩阵中的弹性参数，为横向各向同
性介质的刚度矩阵；
所述横向各向同性介质的刚度矩阵为：
（3）
其中，利用Thomsen参数将刚度矩阵中的弹性参数表示为：
（4）
式中，、、、、、均为弹性参数；、、均为Thomsen参数，为纵波沿横向各向同性介质模型对称轴方向的速度，为横波沿横向各向同性介质
模型对称轴方向的速度；
采用交错网格方法对横向各向同性介质模型进行剖分，并确定正应力、切应力以
及质点振动速度场沿各方向的分量在交错网格上的位置，具体操作如下：将规则节点处设置为一号位置；将在方向位于规则的节点、在y方向上位于
的半节点、在z方向上位于的半节点处设置为二号位置；将在方向位于的
半节点、在y方向上位于规则的的节点、在z方向上位于的半节点处设置为三号位
置；将在方向位于的半节点、在y方向上位于的半节点、在方向位于规则
的节点处设置为四号位置；将在方向位于的半节点、在y方向上位于规则的
的节点、在方向位于规则的节点处设置为五号位置；将在方向位于规则的节点、在
y方向上位于的半节点、在方向位于规则的节点处设置为六号位置；将在方向
位于规则的节点、在y方向上位于规则的的节点、在z方向上位于的半节点处
设置为七号位置。
优选地，所述步骤1 中，将所述正应力、正应力、正应力均设置在一号
位置；切应力设置在二号位置；切应力设置在三号位置；切应力设置在四号位
置；质点振动速度矢量场沿方向的分量设置在五号位置；质点振动速度矢量场沿方
向的分量设置在六号位置；质点振动速度矢量场沿方向的分量设置在七号位置；
对横向各向同性介质弹性波一阶速度-应力方程进行离散化，得到：
（5）
（6）
（7）
（8）
（9）
（10）
（11）
（12）
（13）
式中，为时间采样间隔；为横向方向上的采样间隔，为纵向方向上
的采样间隔，为垂直方向上的采样间隔；为阶交错网格差分系数，为
阶交错网格差分系数的下标，，为半阶数；为整时间网格点，
为半时间网格点。
优选地，所述步骤2中，对横向各向同性介质的刚度矩阵进行分解，将刚度矩阵
分解为纵波刚度矩阵和横波刚度矩阵，如公式（14）和公式（15）所示：
（14）
（15）
式中，为纵波刚度矩阵，为qP波，、、、、、、、、、、、均为纵波弹性参数；为横波刚度矩阵，为qS波，、、、、、、、、、、、均为横波弹性参数；
选取仅与有关的弹性参数作为纵波弹性参数、仅与有关的弹性参数
作为横波弹性参数，则纵波弹性参数和横波弹性参数表示为：
（16）
（17）
此时横波刚度矩阵表示为：
（18）。
优选地，所述步骤3中，引入qS波辅助应力，利用分解得到的横波弹性参数构造
qS波应力，再利用qS波应力计算qS波速度，构建由待定的横波刚度矩阵表示的qS波
解耦一阶速度-应力方程；
所述由待定的横波刚度矩阵表示的qS波解耦一阶速度-应力方程为：
（19）
（20）
式中，为qS波速度，取值为，其中，为qS波速度矢量场沿方向
的分量，为qS波速度矢量场沿方向的分量，为qS波速度矢量场沿方向的分量；为指标索引所对应qS波速度在时间上的一阶导数；为qS波应力，取值为，其中，、、均为正应力，、、均为切
应力；为指标索引所对应的应力对方向求一阶导数；为指标索引所对应的应
力对时间求一阶导数；为公式（18）中横波刚度矩阵的横波弹性参数；
为了消除振幅畸变，对零阶伪亥姆霍兹算子除以纵横波垂直速度的平方差进行振幅修正，如公式（21）所示：
（21）
式中，为修正的零阶伪亥姆霍兹算子；、均为局部各向异性参数；
为横向方向上的偏导数，为纵向方向上的偏导数，为垂直方向上的偏导数；、、均为中间量计算函数；
利用修正的零阶伪亥姆霍兹分解算子，基于横波无散的思想，得到修正后零阶伪
亥姆霍兹算子与qS波速度之间的关系为：
（22）
将公式（19）和公式（20）代入公式（22）中，通过反证法计算得到待定的横波刚度矩阵弹性参数为：
（23）
式中，、、均为计算中间值，其中，取值为，取值为，取值
为；
将公式（23）中计算得到的横波刚度矩阵弹性参数代入公式（20）中，得到横向各向同性介质中qS波解耦的一阶速度-应力方程，如公式（24）所示：
（24）
式中，为时间。
优选地，所述步骤4中，qP波解耦一阶速度-应力方程如公式（25）所示：
（25）
式中，为解耦的qP波应力在横向方向上的分量，为解耦的qP波应力在纵
向方向上的分量，为解耦的qP波应力在垂直方向上的分量；为解耦的qP波质点
振动速度在横向方向上的分量，为解耦的qP波质点振动速度在纵向方向上的分量，为解耦的qP波质点振动速度在垂直方向上的分量。
优选地，所述步骤4中，利用交错网格的有限差分格式求解横向各向同性介质的弹
性波一阶速度-应力方程，计算得到总波场质点振动速度场，，确定总波场
质点振动速度；
将计算得到的总波场质点振动速度代入qP波解耦一阶速度-应力方程中，基于交
错网格的有限差分格式求解解耦qP波质点振动速度，，确定qP波质点振
动速度，得到qP波波场；
将计算得到的总波场质点振动速度代入横向各向同性介质中qS波解耦的一阶速
度-应力方程中，基于交错网格的有限差分格式求解解耦qS波质点振动速度，，确定qS波质点振动速度，得到qS波波场。
本发明所带来的有益技术效果：
本发明提出的一种横向各向同性介质波场分离的横波解耦方程构建方法，针对横向各向同性介质弹性波波场分离中横波波场耦合严重的问题，基于横波无散的反证法在时空域中构建横向各向同性介质的qS波解耦一阶速度-应力方程，结合qP波解耦一阶速度-应力方程，在时空域实现了横向各向同性介质中弹性波场的无串扰分离，该方法实施过程简单，且为后续波场的准确成像奠定了基础。
附图说明
图1为本发明一种横向各向同性介质波场分离的横波解耦方程构建方法的流程图。
图2为本发明横向各向同性介质模型中交错网格的示意图。
图3为x分量波场分离结果的单道波形图。
图4为z分量波场分离结果的单道波形图。
图5为双层介质模型的示意图。
图中：1、第一位置，2、第二位置，3、第三位置，4、第四位置，5、第五位置，6、第六位置，7、第七位置。
具体实施方式
下面结合附图以及具体实施方式对本发明作进一步详细说明：
实施例1
本发明提出了一种横向各向同性介质波场分离的横波解耦方程构建方法，如图1所示，具体包括以下步骤：
步骤1，输入横向各向同性介质模型，对横向各向同性介质弹性波场延拓，具体步骤如下：
将横向各向同性介质弹性波一阶速度-应力方程表示为：
（1）
（2）
式中，为横向各向同性介质模型的横向方向，为横向各向同性介质模型的纵
向方向，为横向各向同性介质模型的垂直方向；为横向各向同性介质的密度；为
指标索引，取值为；代表质点振动速度，取值为,其中，为质点振动
速度矢量场沿方向的分量，为质点振动速度矢量场沿方向的分量，为质点振动速
度矢量场沿方向的分量；为指标索引所对应质点振动速度在时间上的一阶导数；
为指标索引所对应的质点振动速度对方向求一阶导数，为指标索引所对应的质点
振动速度对方向求一阶导数，代表应力，取值为，其
中，、、均为正应力，、、均为切应力；为指标索引所对应的应
力对方向求一阶导数；为指标索引所对应的应力对时间求一阶导数；为横向各
向同性介质刚度矩阵中的弹性参数，为横向各向同性介质的刚度矩阵；
所述横向各向同性介质的刚度矩阵为：
（3）
其中，利用Thomsen参数将刚度矩阵中的弹性参数表示为：
（4）
式中，、、、、、均为弹性参数；、、均为Thomsen参数，为纵波沿横向各向同性介质模型对称轴方向的速度，为横波沿横向各向同性介质
模型对称轴方向的速度。
采用交错网格方法对横向各向同性介质模型进行剖分，建立并确定正应力、切应
力以及质点振动速度场沿各方向的分量在交错网格上的位置，建立如图2所示的交错网格，
具体操作如下：将规则节点处设置为一号位置；将在方向位于规则的节点、
在y方向上位于的半节点、在z方向上位于的半节点处设置为二号位置；将在
方向位于的半节点、在y方向上位于规则的的节点、在z方向上位于的半节点
处设置为三号位置；将在方向位于的半节点、在y方向上位于的半节点、在方
向位于规则的节点处设置为四号位置；将在方向位于的半节点、在y方向上位于
规则的的节点、在方向位于规则的节点处设置为五号位置；将在方向位于规则的节点、在y方向上位于的半节点、在方向位于规则的节点处设置为六号位置；将
在方向位于规则的节点、在y方向上位于规则的的节点、在z方向上位于的半节
点处设置为七号位置。
所述正应力、正应力、正应力均设置在一号位置；切应力设置在二
号位置；切应力设置在三号位置；切应力设置在四号位置；质点振动速度矢量场沿方向的分量设置在五号位置；质点振动速度矢量场沿方向的分量设置在六号位
置；质点振动速度矢量场沿方向的分量设置在七号位置。
对公式（1）和公式（2）所示的横向各向同性介质弹性波一阶速度-应力方程进行离散化，得到：
（5）
（6）
（7）
（8）
（9）
（10）
（11）
（12）
（13）
式中，为时间采样间隔；为横向方向上的采样间隔，为纵向方向上
的采样间隔，为垂直方向上的采样间隔；为阶交错网格差分系数，为
阶交错网格差分系数的下标，，为半阶数；为整时间网格点，为半
时间网格点。
步骤2，基于方程解耦将弹性刚度系数近似分解为待定的纵横波刚度矩阵弹性参数，具体步骤如下：
对横向各向同性介质的刚度矩阵进行分解，将刚度矩阵分解为纵波刚度矩阵和横波刚度矩阵，如公式（14）和公式（15）所示：
（14）
（15）
式中，为纵波刚度矩阵，为qP波，、、、、、、、、、、、均为纵波弹性参数；为横波刚度矩阵，为qS波，、、、、、、、、、、、均为横波弹性参数。
选取仅与有关的弹性参数作为纵波弹性参数、仅与有关的弹性参数
作为横波弹性参数，则纵波弹性参数和横波弹性参数表示为：
（16）
（17）
此时横波刚度矩阵表示为：
（18）。
步骤3，引入修正伪亥姆霍兹算子，基于横波无散的反证法求解待定的横波刚度矩阵弹性参数，利用横波刚度矩阵弹性参数构造qS波应力，根据qS波应力求取qS波质点振动速度，构造qS波解耦一阶速度-应力方程，具体步骤如下：
引入qS波辅助应力，利用分解得到的横波弹性参数构造qS波应力，再利用qS波
应力计算qS波速度，构建由待定的横波刚度矩阵表示的qS波解耦一阶速度-应力方
程；
所述由待定的横波刚度矩阵表示的qS波解耦一阶速度-应力方程为：
（19）
（20）
式中，为qS波速度，取值为，其中，为qS波速度矢量场沿方
向的分量，为qS波速度矢量场沿方向的分量，为qS波速度矢量场沿方向的分量；为指标索引所对应qS波速度在时间上的一阶导数；为qS波应力，取值为，其中，、、均为正应力，、、均为切
应力；为指标索引所对应的应力对方向求一阶导数；为指标索引所对应的应
力对时间求一阶导数；为公式（18）中横波刚度矩阵的横波弹性参数。
为了消除振幅畸变，对零阶伪亥姆霍兹算子除以纵横波垂直速度的平方差进行振幅修正，如公式（21）所示：
（21）
式中，为修正的零阶伪亥姆霍兹算子；、、、、均
为局部各向异性参数；为横向方向上的偏导数，为纵向方向上的偏导数，为垂
直方向上的偏导数；、、均为中间量计算函数。
利用修正的零阶伪亥姆霍兹分解算子，基于横波无散的思想，得到修正后零阶伪
亥姆霍兹算子与qS波速度之间的关系为：
（22）
将公式（19）和公式（20）代入公式（22）中，通过反证法计算得到待定的横波刚度矩阵弹性参数为：
（23）
式中，、、均为计算中间值，其中，取值为，取值为，取值
为。
将公式（23）中计算得到的横波刚度矩阵弹性参数代入公式（20）中，得到横向各向同性介质中qS波解耦的一阶速度-应力方程，如公式（24）所示：
（24）
式中，为时间。
步骤4，将波场延拓得到的质点振动速度场代入qS波解耦一阶速度-应力方程和qP波解耦一阶速度-应力方程中，解耦分离纵横波波场。
其中，所述qP波解耦一阶速度-应力方程如公式（25）所示：
（25）
式中，为解耦的qP波应力在横向方向上的分量，为解耦的qP波应力在纵
向方向上的分量，为解耦的qP波应力在垂直方向上的分量；为解耦的qP波质点
振动速度在横向方向上的分量，为解耦的qP波质点振动速度在纵向方向上的分量，为解耦的qP波质点振动速度在垂直方向上的分量。
利用交错网格的有限差分格式求解公式（1）和公式（2）所示的横向各向同性介质
的弹性波一阶速度-应力方程，计算得到总波场质点振动速度场，。
将计算得到的总波场质点振动速度代入如公式（25）所示的qP波解耦一阶速度-
应力方程中，基于交错网格的有限差分格式求解解耦qP波质点振动速度，，确定qP波质点振动速度，得到qP波波场；
将计算得到的总波场质点振动速度代入如公式（24）所示的横向各向同性介质中
qS波解耦的一阶速度-应力方程中，基于交错网格的有限差分格式求解解耦qS波质点振动
速度，，确定qS波质点振动速度，得到qS波波场。
实施例2
本实施例采用实施例1中所提出的横向各向同性介质波场分离的横波解耦方程构建方法进行波场分离。
本实施例中采用尺寸为的均匀各向异性介质模型，均匀各向异性介质模型的尺寸
设置为4.5km×4.5km，模型密度设置为，模型中的纵波速度设置为，横波
速度设置为，各向异性参数设置为0.15、设置为0.13。将空间采样的横向间隔
和纵向间隔均设置为10m，时间采样间隔设置为1ms，震源采用爆炸震源且设置于模型中心，
震源子波为雷克子波，主频为25Hz。
采用交错网格空间十阶差分，时间二阶差分格式进行数值计算，得到0.5s的波场快照。由波场快照可得，本发明方法能够有效分离出大部分的P波波场和S波波场。
为了进一步验证本发明方法对于波场分离的有效性，抽取x分量(x=1.55km)和z分量（x=1.4km)在网格点处的单道波形结果，得到如图3所示的x分量波场分离结果的单道波形以及如图4所示的z分量波场分离结果的单道波形，由单道波形图能够直观看出残余振幅相对关系（图3和图4中圈中部分均为残余项）。分析图3和图4可得，本实施例中对于x分量，纵波中残余横波的最大相对振幅为1.36%，横波中残余纵波的最大相对振幅为2.11%；对于z分量，纵波中残余横波的最大相对振幅为0.09%，横波中残余纵波的最大相对振幅为0.98%。
由此可见，波场残余的最大振幅与原波场振幅相差了一个数量级以上，从而验证了本发明方法能够有效分离纵波波场和横波波场，且分离后的纵波波场和横波波场之间无串扰。
实施例3
为了验证本发明所提出横向各向同性介质波场分离的横波解耦方程构建方法在界面处的波场分离效果，将实施例1中的横向各向同性介质波场分离的横波解耦方程构建方法应用于双层模型中。
本实施例中双层介质模型如图5所示，模型尺寸设置为4.5km×4.5km，双层介质模
型内设置有两层水平层，两个水平层之间的分界面位于垂向方向2.25km处，其中，第一水平
层的纵波速度设置为、横波速度设置为，第一水平层的密度设置为，各向异性参数设置为0、设置为0；第二水平层的纵波速度设置为、横
波速度设置为，第二水平层的密度设置为，各向异性参数设置为0.2、
设置为0.15。将空间采样的横向间隔和纵向间隔均设置为10m，时间采样间隔设置为1ms，震
源采用爆炸震源，震源子波为雷克子波，主频为25Hz，震源设置于(2.25km，1.85km)处，设置
接收时间为0.55s。
采用交错网格空间十阶差分，时间二阶差分格式进行数值计算得到波场快照，由波场快照中能够明显看出第一水平层与第二水平层分界面处分离出的纵波波场与横波波场之间无串扰，验证了本发明方法能够实现对相邻介质分界面处纵波波场和横波波场的无串扰分离。
实施例4
本实施例中采用横向各向同性介质模型验证本发明方法的准确性。
本实施例中横向各向同性介质模型的横向宽度设置为40km、纵向深度设置为6.7km，横向网格间距设置为25m、纵向网格间距设置为10m，横向各向同性介质模型表面均匀设置有800个检波点和800个炮点，相邻检波点和相邻炮点之间的间隔设置50m，震源子波为Ricker子波，主频为20Hz。
利用本发明方法获取横向各向同性介质模型的纵波速度、横波速度以及各向异性参数模型，得到四种弹性波成像结果。分析弹性波成像结果可得，四种成像结果均可准确成像且无明显串扰假象，成像深度与横向各向同性介质模型实际深度一致，且能够准确获取断层、小构造、高陡等构造的成像位置和形态，即便在深层也能得到清晰的成像结果，具有较高的分辨率以及信噪比，从而验证了本发明方法的准确性。
当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护范围。
Claims (3)
1.一种横向各向同性介质波场分离的横波解耦方程构建方法，其特征在于，具体包括以下步骤：
步骤1，输入横向各向同性介质模型，对横向各向同性介质进行弹性波场延拓；
步骤2，基于方程解耦将弹性刚度系数近似分解为待定的纵横波刚度矩阵弹性参数；
步骤3，引入修正伪亥姆霍兹算子，基于横波无散的反证法求解待定的横波刚度矩阵弹性参数，利用横波刚度矩阵弹性参数构造qS波应力，根据qS波应力求取qS波质点振动速度，构造qS波解耦一阶速度-应力方程；
步骤4，将波场延拓得到的质点振动速度场代入qS波解耦一阶速度-应力方程和qP波解耦一阶速度-应力方程中，解耦分离纵横波波场；
所述步骤1中，将横向各向同性介质弹性波一阶速度-应力方程表示为：
（1）
（2）
式中，为横向各向同性介质模型的横向方向，/>为横向各向同性介质模型的纵向方向，为横向各向同性介质模型的垂直方向；/>为横向各向同性介质的密度；/>为指标索引，取值为/>、/>、/>；/>代表质点振动速度，取值为/>、/>、/>，其中，/>为质点振动速度矢量场沿/>方向的分量，/>为质点振动速度矢量场沿/>方向的分量，/>为质点振动速度矢量场沿/>方向的分量；/>为指标索引/>所对应质点振动速度在时间上的一阶导数；/>为指标索引/>所对应的质点振动速度对/>方向求一阶导数，/>为指标索引/>所对应的质点振动速度对/>方向求一阶导数，/>代表应力，取值为/>、/>、/>、/>、/>、/>，其中，/>、/>、/>均为正应力，、/>、/>均为切应力；/>为指标索引/>所对应的应力对/>方向求一阶导数；/>为指标索引/>所对应的应力对时间求一阶导数；/>为横向各向同性介质刚度矩阵中的弹性参数，/>为横向各向同性介质的刚度矩阵；
所述横向各向同性介质的刚度矩阵为：
（3）
其中，利用Thomsen参数将刚度矩阵中的弹性参数表示为：
（4）
式中，、/>、/>、/>、/>、/>均为弹性参数；/>、/>、/>均为Thomsen参数，/>为纵波沿横向各向同性介质模型对称轴方向的速度，/>为横波沿横向各向同性介质模型对称轴方向的速度；
采用交错网格方法对横向各向同性介质模型进行剖分，并确定正应力、切应力以及质点振动速度场沿各方向的分量在交错网格上的位置，具体操作如下：将规则节点处设置为一号位置；将在/>方向位于规则的/>节点、在y方向上位于/>的半节点、在z方向上位于/>的半节点处设置为二号位置；将在/>方向位于/>的半节点、在y方向上位于规则的/>的节点、在z方向上位于/>的半节点处设置为三号位置；将在/>方向位于/>的半节点、在y方向上位于/>的半节点、在/>方向位于规则的/>节点处设置为四号位置；将在/>方向位于/>的半节点、在y方向上位于规则的/>的节点、在/>方向位于规则的/>节点处设置为五号位置；将在/>方向位于规则的/>节点、在y方向上位于/>的半节点、在/>方向位于规则的/>节点处设置为六号位置；将在/>方向位于规则的/>节点、在y方向上位于规则的/>的节点、在z方向上位于/>的半节点处设置为七号位置；
所述步骤1 中，将所述正应力、正应力/>、正应力/>均设置在一号位置；切应力/>设置在二号位置；切应力/>设置在三号位置；切应力/>设置在四号位置；质点振动速度矢量场沿/>方向的分量/>设置在五号位置；质点振动速度矢量场沿/>方向的分量/>设置在六号位置；质点振动速度矢量场沿/>方向的分量/>设置在七号位置；
对横向各向同性介质弹性波一阶速度-应力方程进行离散化，得到：
（5）
（6）
（7）
（8）
（9）
（10）
（11）
（12）
（13）
式中，为时间采样间隔；/>为横向方向/>上的采样间隔，/>为纵向方向/>上的采样间隔，/>为垂直方向/>上的采样间隔；/>为/>阶交错网格差分系数，/>为/>阶交错网格差分系数的下标，/>，/>为半阶数；/>为整时间网格点，/>为半时间网格点；
所述步骤2中，对横向各向同性介质的刚度矩阵进行分解，将刚度矩阵/>分解为纵波刚度矩阵/>和横波刚度矩阵/>，如公式（14）和公式（15）所示：
（14）
（15）
式中，为纵波刚度矩阵，/>为qP波，/>、/>、/>、/>、/>、/>、/>、/>、/>、/>、/>、均为纵波弹性参数；/>为横波刚度矩阵，/>为qS波，/>、/>、/>、/>、/>、/>、/>、/>、、/>、/>、/>均为横波弹性参数；
选取仅与有关的弹性参数作为纵波弹性参数、仅与/>有关的弹性参数作为横波弹性参数，则纵波弹性参数和横波弹性参数表示为：
（16）
（17）
此时横波刚度矩阵表示为：
（18）；
所述步骤3中，引入qS波辅助应力，利用分解得到的横波弹性参数构造qS波应力，再利用qS波应力计算qS波速度/>，构建由待定的横波刚度矩阵/>表示的qS波解耦一阶速度-应力方程；
所述由待定的横波刚度矩阵表示的qS波解耦一阶速度-应力方程为：
（19）
（20）
式中，为qS波速度，取值为/>、/>、/>，其中，/>为qS波速度矢量场沿/>方向的分量，为qS波速度矢量场沿/>方向的分量，/>为qS波速度矢量场沿/>方向的分量；/>为指标索引/>所对应qS波速度在时间上的一阶导数；/>为qS波应力，取值为/>、/>、/>、/>、/>、，其中，/>、/>、/>均为正应力，/>、/>、/>均为切应力；/>为指标索引/>所对应的应力对/>方向求一阶导数；/>为指标索引/>所对应的应力对时间求一阶导数；/>为公式（18）中横波刚度矩阵/>的横波弹性参数；
为了消除振幅畸变，对零阶伪亥姆霍兹算子除以纵横波垂直速度的平方差进行振幅修正，如公式（21）所示：
（21）
式中，为修正的零阶伪亥姆霍兹算子；/>、/>均为局部各向异性参数；/>为横向方向/>上的偏导数，/>为纵向方向/>上的偏导数，/>为垂直方向/>上的偏导数；/>、/>、均为中间量计算函数；
利用修正的零阶伪亥姆霍兹分解算子，基于横波无散的思想，得到修正后零阶伪亥姆霍兹算子与qS波速度之间的关系为：
（22）
将公式（19）和公式（20）代入公式（22）中，通过反证法计算得到待定的横波刚度矩阵弹性参数为：
（23）
式中，、/>、/>均为计算中间值，其中，/>取值为/>，/>取值为/>，/>取值为/>；
将公式（23）中计算得到的横波刚度矩阵弹性参数代入公式（20）中，得到横向各向同性介质中qS波解耦的一阶速度-应力方程，如公式（24）所示：
（24）
式中，为时间。
2.根据权利要求1所述的横向各向同性介质波场分离的横波解耦方程构建方法，其特征在于，所述步骤4中，qP波解耦一阶速度-应力方程如公式（25）所示：
（25）
式中，为解耦的qP波应力在横向方向/>上的分量，/>为解耦的qP波应力在纵向方向上的分量，/>为解耦的qP波应力在垂直方向/>上的分量；/>为解耦的qP波质点振动速度在横向方向/>上的分量，/>为解耦的qP波质点振动速度在纵向方向/>上的分量，/>为解耦的qP波质点振动速度在垂直方向/>上的分量。
3.根据权利要求2所述的横向各向同性介质波场分离的横波解耦方程构建方法，其特征在于，所述步骤4中，利用交错网格的有限差分格式求解横向各向同性介质的弹性波一阶速度-应力方程，计算得到总波场质点振动速度场，/>，确定总波场质点振动速度；
将计算得到的总波场质点振动速度代入qP波解耦一阶速度-应力方程中，基于交错网格的有限差分格式求解解耦qP波质点振动速度/>，/>，确定qP波质点振动速度，得到qP波波场；
将计算得到的总波场质点振动速度代入横向各向同性介质中qS波解耦的一阶速度-应力方程中，基于交错网格的有限差分格式求解解耦qS波质点振动速度/>，/>，确定qS波质点振动速度，得到qS波波场。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311316458.3A CN117075197B (zh) | 2023-10-12 | 2023-10-12 | 一种横向各向同性介质波场分离的横波解耦方程构建方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311316458.3A CN117075197B (zh) | 2023-10-12 | 2023-10-12 | 一种横向各向同性介质波场分离的横波解耦方程构建方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN117075197A CN117075197A (zh) | 2023-11-17 |
| CN117075197B true CN117075197B (zh) | 2024-02-06 |
Family
ID=88711911
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202311316458.3A Active CN117075197B (zh) | 2023-10-12 | 2023-10-12 | 一种横向各向同性介质波场分离的横波解耦方程构建方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN117075197B (zh) |
Citations (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2012036689A1 (en) * | 2010-09-16 | 2012-03-22 | Halliburton Energy Services, Inc. | Combined sonic/pulsed neutron cased hole logging tool |
| CN104133241A (zh) * | 2014-07-31 | 2014-11-05 | 中国科学院地质与地球物理研究所 | 波场分离方法和装置 |
| WO2016083893A1 (en) * | 2014-11-25 | 2016-06-02 | Cgg Services Sa | Modeling an elastic stiffness tensor in a transverse isotropic subsurface medium |
| CN112904426A (zh) * | 2021-03-27 | 2021-06-04 | 中国石油大学(华东) | 一种解耦弹性波逆时偏移方法、系统及应用 |
| CN113406698A (zh) * | 2021-05-24 | 2021-09-17 | 中国石油大学(华东) | 一种基于纵横波解耦的双相介质弹性波逆时偏移成像方法 |
| CN113885079A (zh) * | 2021-08-23 | 2022-01-04 | 中国石油大学(华东) | 基于弹性波场解耦的高精度多方位逆时震源成像方法 |
| CN115657120A (zh) * | 2022-10-21 | 2023-01-31 | 中国石油大学(华东) | 一种弹性波应力张量双点积地震成像方法、装置及设备 |
| CN115685334A (zh) * | 2022-10-24 | 2023-02-03 | 中国石油大学(华东) | 一种各向异性弹性波场分解方法、装置及计算机设备 |
| CN116068621A (zh) * | 2021-11-01 | 2023-05-05 | 中国石油天然气股份有限公司 | 一种基于刚度矩阵分解的各向异性介质正演方法及系统 |
| CN116520418A (zh) * | 2023-04-13 | 2023-08-01 | 中国石油大学(华东) | 一种弹性波角度域共成像点道集高效提取方法 |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| AU2011224165B2 (en) * | 2010-03-12 | 2013-10-10 | Cggveritas Services (Us) Inc. | Methods and systems for performing azimuthal simultaneous elastic inversion |
| US11360224B2 (en) * | 2019-05-03 | 2022-06-14 | Exxonmobil Upstream Research Company | Inversion, migration, and imaging related to isotropic wave-mode- independent attenuation |
-
2023
- 2023-10-12 CN CN202311316458.3A patent/CN117075197B/zh active Active
Patent Citations (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2012036689A1 (en) * | 2010-09-16 | 2012-03-22 | Halliburton Energy Services, Inc. | Combined sonic/pulsed neutron cased hole logging tool |
| CN104133241A (zh) * | 2014-07-31 | 2014-11-05 | 中国科学院地质与地球物理研究所 | 波场分离方法和装置 |
| WO2016083893A1 (en) * | 2014-11-25 | 2016-06-02 | Cgg Services Sa | Modeling an elastic stiffness tensor in a transverse isotropic subsurface medium |
| CN112904426A (zh) * | 2021-03-27 | 2021-06-04 | 中国石油大学(华东) | 一种解耦弹性波逆时偏移方法、系统及应用 |
| CN113406698A (zh) * | 2021-05-24 | 2021-09-17 | 中国石油大学(华东) | 一种基于纵横波解耦的双相介质弹性波逆时偏移成像方法 |
| CN113885079A (zh) * | 2021-08-23 | 2022-01-04 | 中国石油大学(华东) | 基于弹性波场解耦的高精度多方位逆时震源成像方法 |
| CN116068621A (zh) * | 2021-11-01 | 2023-05-05 | 中国石油天然气股份有限公司 | 一种基于刚度矩阵分解的各向异性介质正演方法及系统 |
| CN115657120A (zh) * | 2022-10-21 | 2023-01-31 | 中国石油大学(华东) | 一种弹性波应力张量双点积地震成像方法、装置及设备 |
| CN115685334A (zh) * | 2022-10-24 | 2023-02-03 | 中国石油大学(华东) | 一种各向异性弹性波场分解方法、装置及计算机设备 |
| CN116520418A (zh) * | 2023-04-13 | 2023-08-01 | 中国石油大学(华东) | 一种弹性波角度域共成像点道集高效提取方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN117075197A (zh) | 2023-11-17 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN112904426B (zh) | 一种解耦弹性波逆时偏移方法、系统及应用 | |
| CN101201409B (zh) | 一种地震数据变相位校正方法 | |
| CN111158049B (zh) | 一种基于散射积分法的地震逆时偏移成像方法 | |
| CN109765616B (zh) | 一种保幅波场延拓校正方法及系统 | |
| MX2011003850A (es) | Estimado de señal de dominio de imagen a interferencia. | |
| CN113885079B (zh) | 基于弹性波场解耦的高精度多方位逆时震源成像方法 | |
| CN106019376A (zh) | 一种频率驱动空变q值模型构建的地震波补偿方法 | |
| Ning et al. | Numerical investigation of Rayleigh-wave propagation on canyon topography using finite-difference method | |
| CN104570095A (zh) | 一种基于Radon变换消除斜缆虚反射的方法 | |
| CN108919356A (zh) | 一种稳定的大沙漠区衰减补偿逆时偏移成像系统及方法 | |
| CN117075197B (zh) | 一种横向各向同性介质波场分离的横波解耦方程构建方法 | |
| CN103558636A (zh) | 一种从叠后地震数据采集脚印衰减的方法 | |
| Foti et al. | Spatial sampling issues in fk analysis of surface waves | |
| CN105182414A (zh) | 一种基于波动方程正演去除直达波的方法 | |
| CN107229066A (zh) | 基于地面地震构造约束的vsp数据全波形反演建模方法 | |
| CN110780346A (zh) | 一种隧道超前探测复杂地震波场的分离方法 | |
| CN113536638B (zh) | 基于间断有限元和交错网格的高精度地震波场模拟方法 | |
| Gao et al. | Acquisition and processing pitfall with clipped traces in surface-wave analysis | |
| CN116068621A (zh) | 一种基于刚度矩阵分解的各向异性介质正演方法及系统 | |
| CN113139266B (zh) | 纵、横波数值模拟方法及系统 | |
| MX2011003852A (es) | Atributos de procesamiento de imagen con inversion de tiempo. | |
| Yang et al. | Imaging dispersion of MASW data with the rammer source | |
| CN116413801B (zh) | 一种各向异性介质弹性波高精度成像方法和系统 | |
| CN115453621B (zh) | 基于一阶速度-应力方程的纵横波解耦分离假象去除方法 | |
| Huang et al. | Observed evolution of linear and nonlinear effects at the Dahan downhole array, Taiwan: Analysis of the September 21, 1999 M 7.3 Chi-Chi earthquake sequence |
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 |