CN116755148A - 正交各向异性介质多方位反射波走时反演方法 - Google Patents
正交各向异性介质多方位反射波走时反演方法 Download PDFInfo
- Publication number
- CN116755148A CN116755148A CN202310595476.3A CN202310595476A CN116755148A CN 116755148 A CN116755148 A CN 116755148A CN 202310595476 A CN202310595476 A CN 202310595476A CN 116755148 A CN116755148 A CN 116755148A
- Authority
- CN
- China
- Prior art keywords
- coordinate system
- inversion
- medium
- travel time
- azimuth
- 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
Links
- 238000000034 method Methods 0.000 title claims abstract description 69
- 238000006243 chemical reaction Methods 0.000 claims abstract description 17
- 238000011156 evaluation Methods 0.000 claims abstract description 4
- 239000011159 matrix material Substances 0.000 claims description 97
- 239000013598 vector Substances 0.000 claims description 27
- 230000008569 process Effects 0.000 claims description 14
- 230000008859 change Effects 0.000 claims description 10
- 238000004364 calculation method Methods 0.000 claims description 9
- 238000010586 diagram Methods 0.000 claims description 6
- 230000009466 transformation Effects 0.000 claims description 6
- 230000015572 biosynthetic process Effects 0.000 claims description 5
- 230000003044 adaptive effect Effects 0.000 claims description 3
- 238000004458 analytical method Methods 0.000 claims description 3
- 238000007619 statistical method Methods 0.000 claims description 3
- 238000011144 upstream manufacturing Methods 0.000 claims description 3
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 claims 1
- 230000000644 propagated effect Effects 0.000 claims 1
- 239000004576 sand Substances 0.000 claims 1
- 208000037516 chromosome inversion disease Diseases 0.000 description 76
- 238000003384 imaging method Methods 0.000 description 5
- 239000012530 fluid Substances 0.000 description 4
- 238000012545 processing Methods 0.000 description 4
- 238000012795 verification Methods 0.000 description 4
- 239000011435 rock Substances 0.000 description 3
- 238000010276 construction Methods 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- 238000013508 migration Methods 0.000 description 2
- 230000005012 migration Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 239000003208 petroleum Substances 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 239000003079 shale oil Substances 0.000 description 2
- 230000008901 benefit Effects 0.000 description 1
- 238000004422 calculation algorithm Methods 0.000 description 1
- 238000012512 characterization method Methods 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 239000013078 crystal Substances 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 235000011850 desserts Nutrition 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000005286 illumination Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 229910052500 inorganic mineral Inorganic materials 0.000 description 1
- 239000011707 mineral Substances 0.000 description 1
- 230000002265 prevention Effects 0.000 description 1
- 229920006395 saturated elastomer Polymers 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. analysis, for interpretation, for correction
- 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
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/16—Receiving elements for seismic signals; Arrangements or adaptations of receiving elements
- G01V1/18—Receiving elements, e.g. seismometer, geophone or torque detectors, for localised single point measurements
- G01V1/181—Geophones
-
- 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. analysis, for interpretation, for correction
- G01V1/282—Application of seismic models, synthetic seismograms
Abstract
本发明提供一种正交各向异性介质多方位反射波走时反演方法，包括：步骤1、建立全局坐标系及局部坐标系；步骤2、确立局部坐标系下正交各向异性介质单方位反射P波走时近似公式；步骤3、进行基于单方位走时的局部坐标系正交各向异性介质各向异性参数反演；步骤4、建立全局坐标系下的各向异性参数转换公式；步骤5、基于转换公式进行全局各向异性参数反演；步骤6、进行反演结果精度评价。该正交各向异性介质多方位反射波走时反演方法能够更加充分的利用多方位地震数据，提供更加丰富的各向异性信息，为后续的地层属性识别及地质解释工作奠定了基础。
Description
技术领域
本发明涉及石油地球物理勘探技术领域，特别是涉及到一种正交各向异性介质多方位反射波走时反演方法。
背景技术
地下介质中广泛存在各向异性，在地震勘探中，地下介质的各向异性主要表现为速度各向异性。地下介质的各向异性成因多种多样，包括晶体矿物的定向排列，沉积岩精细分层，裂缝定向排列及地应力影响等。在页岩油气勘探开发中，裂缝控制油气运移的通道与储集空间，是地质“甜点”的重要识别标志；在水平压裂过程中，为保证压裂效果，需要对地层水平最大主应力方向以及水平应力差进行识别。因此，在页岩油气勘探领域，方位各向异性的研究越来越受到重视。
具有方位各向异性特征的各向异性介质，根据对称性由简单到复杂可以分为水平横向各向同性(Horizontal Transverse Isotropic,HTI)介质、正交各向异性介质和单斜各向异性介质等。其中，当地下介质中存在两组相互正交的垂直裂隙时，等效的介质为正交各向异性介质。与TI介质相比，正交各向异性介质需要更多各向异性参数进行描述，在TI介质中，由于其存在一个对称轴，因此根据对称轴确定坐标系后，仅需确定其在沿对称轴方向及沿与对称轴方向垂直方向的速度及其在与对称轴平行的轴面内速度的变化梯度。而对于正交各向异性介质，其存在两个对称面，因此根据对称面确定坐标系后，需要确定其沿三个坐标轴方向传播的速度以及其在三个轴面内速度的变化梯度。
关于正交各向异性介质的研究于1997年开始，Tsvakin提出了正交各向异性介质理论并基于其弹性系数矩阵推导了正交各向异性介质的P波速度及相关的各向异性参数。后续关于正交各向异性介质的研究主要包括对反射波走时及的振幅的正演和反演等。走时反演方面主要以NMO速度反演为主，NMO速度是对反射双曲线走时进行正常时差校正时使用的速度。NMO速度随方位的变化较为规律，呈椭圆属性，对NMO速度的反演可以较精确地确定各向异性介质的方位属性。但NMO速度的反演流程更加依赖于先验信息如反射层厚度，垂直方向传播的速度等，因此具有一定的局限性。
在申请号：CN202010817813.5的中国专利申请中，涉及到一种正交各向异性速度反演的方法及系统，方法包括：采集成像点道集和叠加数据；利用模型道与成像点道集中的每一地震数据道在用户定义的处理时窗内做互相关计算；对当前地震数据的每一个互相关时窗做互相关计算；对当前成像点道集中的所有地震数据道进行处理；根据当前成像点道集的剩余时差数据和相关系数数据及叠前时间偏移速度场数据，利用加权最小平方算法计算得到当前成像点道集在所有零炮检距数据样点上的HTI各向异性速度参数；根据HTI各向异性速度及剩余时差的残差部分得到VTI各向异性速度参数；对每个成像点道集进行处理，得到所有地震数据的正交各向异性速度反演结果，该方法具有高效的计算效率，且计算结果精度高。
在申请号：CN201210371061.X的中国专利申请中，涉及到一种地震各向异性参数全波形反演方法及装置，属于石油地球物理勘探中地震各向异性参数预测领域。所述方法包括以下步骤：(1)获取将要进行各向异性参数反演的地震数据，即观测波场数据；(2)对步骤(1)得到的地震数据进行去噪处理得到去噪后的数据；(3)使用步骤(2)得到的去噪后的数据，根据炮点和检波点的坐标，抽取共中心点道集得到地震CMP道集资料，然后利用所述地震CMP道集资料计算出水平地层的层速度；(4)利用步骤(3)得到的水平地层的层速度，通过插值的方式构建用于反演的初始模型；(5)对初始模型的各参数进行微小扰动，生成参数扰动后的模型。
在申请号：CN202011148799.0的中国专利申请中，涉及到一种正交各向异性参数计算方法、装置、电子设备及介质。正交各向异性参数计算方法包括：提取方位道集的剩余时差；利用所述剩余时差反演HTI介质的各向异性参数；对所述方位道集进行时差校正；利用振幅反演VTI介质的各向异性参数。该发明从方位叠前道集出发，将HTI介质与VTI介质分别考虑，分别反演，解决了常规的一次反演中相互干扰，相互影响，进而造成反演结果不准确的问题，为稳定、快速且具有高信噪比的开展两组裂缝各向异性预测技术开辟了另外一条道路。
在申请号：CN201911059377.3的中国专利申请中，涉及到一种正交各向异性介质流体因子与裂缝参数反演方法。包括以下步骤：使用干岩石背景下垂直横向各向同性介质，并根据裂缝弱度推导出正交对称各向异性裂隙介质的刚度矩阵；利用弱各向异性和小裂缝弱度的假设，结合各向异性Gassmann方程导出正交介质中饱和流体岩石弱各向异性近似刚度的新表达式；结合刚度扰动和散射理论，得到正交对称弱各向异性介质中流体和裂隙参数解耦的PP波线性化反射系数；以测井信息作为先验信息，在贝叶斯框架下利用部分入射角叠加方位地震数据实现随偏移距和方位角变化的弹性阻抗叠前反演。该发明能够流体识别和裂隙表征提供可靠的结果，为地震波的传播研究、油气开发和地震灾害预防提供了技术支持。
以上现有技术均与本发明有较大区别，未能解决我们想要解决的技术问题，为此我们发明了一种新的正交各向异性介质多方位反射波走时反演方法。
发明内容
本发明的目的是提供一种能够更加充分的利用多方位地震数据，提供更加丰富的各向异性信息的正交各向异性介质多方位反射波走时反演方法。
本发明的目的可通过如下技术措施来实现：正交各向异性介质多方位反射波走时反演方法，该正交各向异性介质多方位反射波走时反演方法包括：
步骤1、建立全局坐标系及局部坐标系；
步骤2、确立局部坐标系下正交各向异性介质单方位反射P波走时近似公式；
步骤3、进行基于单方位走时的局部坐标系正交各向异性介质各向异性参数反演；
步骤4、建立全局坐标系下的各向异性参数转换公式；
步骤5、基于转换公式进行全局各向异性参数反演；
步骤6、进行反演结果精度评价。
本发明的目的还可通过如下技术措施来实现：
在步骤1，假设一笛卡尔坐标系作为全局坐标系，其x3轴是垂直方向的，x1，x2轴是水平方向且相互正交，轴面(x1，x2)、(x1，x3)分别与介质的两个对称面平行；背景介质为一均匀的地层，地层厚度为H，地层底部的水平反射面与(x1，x2)平面平行；在这种情况下，P波的入射下行和反射上行射线位于垂直于反射面的平面内，并且相对于反射面法线对称。
在步骤1，针对地震波观测系统，定义局部坐标系，反射射线和两个使用的坐标系示意图：全局坐标系xi和局部坐标系x′i；角是震源检波器剖面的方位角，SO是入射射线，OR是反射射线。
在步骤1，地震波从全局坐标系的原点处的震源S激发，在点O处反射，传播至检波器R；将点S和R所在的轴定义为局部坐标系的x′1轴，与x′1轴在同一水平面且与x′1轴垂直的轴为x′2轴，垂直于x′1轴、x′2轴的轴为x′3轴；局部坐标系的原点与震源S一致，并且还表示全局坐标系xi的原点，垂直轴x3在局部和全局两个坐标系中也是相同的；反射射线位于局部坐标系的(x′1，x′3)平面内，与全局坐标系中的坐标平面(x1，x3)夹角为
在步骤2，在二维坐标系下，反射波走时T具有以下形式：
其中，x为偏移距即震源S和检波器R之间的距离，H是水平反射面的深度，与地层的对称平面重合；T＝T(x)表示反射P波的走时；v＝v(n)表示群速度，n为相位矢量；利用时差分析中广泛应用的表示方法：
其中，为归一化的偏移距，α为参考各向同性介质中的P波相速度，T0为参考零偏移距双程走时；公式(1)可改写为：
在步骤2，采用相速度平方近似，即v2(N)～c2(N)近似表示群速度v(N)；在二维坐标系下，正交各向异性介质的群速度可使用各向异性参数表示为：
公式(3)中，N为速度方向矢量，可以通过表示，
涉及的各向异性参数包括εx，εz，ηy，其中，εx，εz分别代表介质在坐标轴x1、x3方向P波相速度与参考速度α的差异；ηy代表介质在坐标平面(x1,x3)内P波速度的变化梯度；局部坐标系下的各向异性参数与局部坐标系下的刚度系数有关：
其中，Aij为地层介质在局部坐标系下的密度归一化刚度系数矩阵中的元素，α为参考各向同性介质中的P波相速度；
经简化，公式(3)可写为以下形式：
步骤3包括：
步骤31，基于零偏移距走时计算εz；
步骤32，构建反演行列式；
步骤33，求解反演行列式。
在步骤31，在零偏移距处，即x＝0时，公式(7)和公式(8)变为：
因此，可通过公式(12)直接求出A参数εz：
在步骤32，根据公式(10)，公式(7)和公式(8)可改写为：
根据公式(11)构建反演行列式如下：
Gm＝d, (12)
其中，G为一N×M的系数矩阵，N为检波器的数量，M为涉及的A参数的数量，此节中M＝2；矩阵G中的行向量为以下形式：
m由2个参数组成，其形式如下：
m＝(εx，ηy)T， (14)
向量d的形式如下：
为归一化偏移距，Tobs为检波器接收到的走时。
在步骤33，公式(12)可以通过最小二乘的方法进行求解，反演的目标函数为：
||Gm-d||2＜κ， (16)
其中，κ为控制反演精度的参数；通过公式(17)对公式(16)进行求解：
其中，为反演算子
其中，GT为矩阵G的转置，I为N×N的单位矩阵，ε为稳定因子，ε越大，反演流程越稳定，但精度较低。
在步骤4，在三维坐标系下，首先假设在平面上沿不同方向分布的地震测线，测线的走向角为由以测线走向为x1轴方向的局部坐标系向全局坐标系转换，转换矩阵R为：
将局部坐标系中的正交各向异性介质各向异性参数使用上标′来与全局坐标系中的正交各向异性介质各向异性参数进行区分；根据Bond变换中的6×6刚度系数矩阵转换法则，局部坐标系中的正交各向异性介质各向异性参数ε′x，η′y和ε′z可以通过以下公式转换为全局坐标系：
ε′z＝εz, (22)
局部坐标系下的各向异性参数定义见公式(6)，全局各向异性参数包括εx，εy，εz，ηx，ηy，ηz，其中，εx，εy，εz分别代表介质在坐标轴x1、x2、x3方向P波相速度与参考速度α的差异；ηx，ηy，ηz分别代表介质在坐标平面(x2,x3)、(x1,x3)、(x1,x2)内P波速度的变化梯度；全局坐标系下的各向异性参数与全局坐标系下的刚度系数有关：
其中，Aij为地层介质在全局坐标系下的密度归一化刚度系数矩阵中的元素，α为参考各向同性介质中的P波相速度。
步骤5包括：
步骤51，构建反演行列式；
步骤52，进行反演行列式求解。
在步骤51，根据公式(20)～(22)，由于沿不同方向分布的地震测线所在局部坐标系之间仅存在方位角的差异，因此各局部坐标系下的A参数εz与全局坐标系下相同，不需要进一步反演；对于剩余五个非零A参数，可通过构建以下行列式进行下一步反演；根据公式(20)和(21)，构建反演行列式如下：
G1m1＝d1, (24)
G2m2＝d2, (25)
其中，G1为N×3的系数矩阵，N为沿不同方向分布的地震测线的数量，由于公式(20)中非零A参数个数为3个，因此矩阵G1的列数为3；G2为N×2的系数矩阵，行数N为沿不同方向分布的地震测线的数量，列数为公式(21)中非零A参数的个数；矩阵G1、G2中的行向量分别为以下形式：
向量m1、m2的形式分别如下：
m1＝(εx，εy，ηz)T， (28)
m2＝(ηx，ηy)T， (29)
向量d1、d2的形式分别如下：
为了使公式(24)、(25)的反问题为适定反问题，需要系数矩阵G1和G2为正定或者超定矩阵；因此，通过多个方位地震测线所在局部坐标系下的A参数反演得到全局坐标系下的A参数，需要至少3个方位的地震反射波走时数据。
在步骤52，公式(24)(25)可以通过最小二乘的方法进行求解，反演的目标函数为：
||G1m1-d1||2＜κ， (32)
||G2m2-d2||2＜κ， (33)
其中，κ为控制反演精度的参数；通过公式(34)对公式(32)(33)进行求解：
其中，为反演算子
其中，GT为矩阵G的转置，I为N×N的单位矩阵，ε为稳定因子，ε越大，反演流程越稳定，但精度较低。
在步骤6，通过将数据协方差矩阵Gd转化为模型协方差矩阵Cm的方法评估所求各向异性参数的精度：
Cm＝HCdHT， (36)
此处模型协方差矩阵Cm表示所求的模型参数向量m的随机分布概率，若随机性不存在，则协方差矩阵为0，此时模型参数向量m没有误差。HT为矩阵H的转置。
在步骤6，如果准确的模型协方差矩阵通过以下方法进行计算：
Cd～σ2I， (37)
其中I为N×N的单位矩阵，σ通过χ2残差统计方法进行计算:
其中： 为公式(17)的计算结果，v为矩阵G行数与列数的差。
本发明中的正交各向异性介质多方位反射波走时反演方法，针对多方位的反射P波走时，首先通过单方位的反射波走时反演单方位的正交各向异性介质各向异性参数，随后根据不同方位的各向异性参数，反演正交各向异性介质在全局坐标系下的各向异性参数的方法。
本发明的优点是提供了一种更可靠、提供更多信息的正交各向异性介质反演方法。在反演过程中减少了对地层厚度及垂直传播速度等先验信息的依赖，同时通过反演能够更多的各向异性参数。从而能够更加充分的利用多方位地震数据，提供更加丰富的各向异性信息，为后续的地层属性识别及地质解释工作奠定了基础。地球物理数据处理人员可以根据本方法处理多方位反射波走时，或依据本发明提供的方法进一步开发PSV波的反演方法。
附图说明
图1为本发明中定义的反射射线和全局及局部坐标系示意图；
图2为本发明的实施例1中的地震测线的分布范围图；
图3为本发明的实施例1中评估所求各向异性参数精度的协方差矩阵；
图4为本发明的实施例1中反演得到的各向异性参数与精确值的对比图；
图5为本发明的实施例2中的地震测线的分布范围图；
图6为本发明的实施例2中评估所求各向异性参数精度的协方差矩阵；
图7为本发明的实施例2中反演得到的各向异性参数与精确值的对比图；
图8为本发明的实施例3中的地震测线的分布范围图；
图9为本发明的实施例3中评估所求各向异性参数精度的协方差矩阵；
图10为本发明的实施例3中反演得到的各向异性参数与精确值的对比图；
图11为本发明的实施例4中的地震测线的分布范围图；
图12为本发明的实施例4中评估所求各向异性参数精度的协方差矩阵；
图13为本发明的实施例4中反演得到的各向异性参数与精确值的对比图；
图14为本发明正交各向异性介质多方位反射波走时反演方法的流程图。
具体实施方式
应该指出，以下详细说明都是示例性的，旨在对本发明提供进一步的说明。除非另有指明，本文使用的所有技术和科学术语具有与本发明所属技术领域的普通技术人员通常理解的相同含义。
需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本发明的示例性实施方式。如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，还应当理解的是，当在本说明书中使用术语“包含”和/或“包括”时，其指明存在特征、步骤、操作和/或它们的组合。
如图14所示，图14为本发明的正交各向异性介质多方位反射波走时反演方法的流程图。该正交各向异性介质多方位反射波走时反演方法包括：
步骤1、建立全局坐标系及局部坐标系；
首先假设一笛卡尔坐标系作为全局坐标系，其x3轴是垂直方向的，x1，x2轴是水平方向且相互正交，轴面(x1，x2)、(x1，x3)分别与介质的两个对称面平行。背景介质为一均匀的地层，地层厚度为H，地层底部的水平反射面与(x1，x2)平面平行。在这种情况下，P波的入射下行和反射上行射线位于垂直于反射面的平面内，并且相对于反射面法线对称。
针对地震波观测系统，定义局部坐标系。如图1所示，反射射线和两个使用的坐标系示意图：全局坐标系xi和局部坐标系x′i。角是震源检波器剖面的方位角，SO是入射射线，OR是反射射线。
地震波从全局坐标系的原点处的震源S激发，在点O处反射，传播至检波器R。将点S和R所在的轴定义为局部坐标系的x′1轴，与x′1轴在同一水平面且与x′1轴垂直的轴为x′2轴，垂直于x′1轴、x′2轴的轴为x′3轴。局部坐标系的原点与震源S一致，并且还表示全局坐标系xi的原点，垂直轴x3在局部和全局两个坐标系中也是相同的。反射射线位于局部坐标系的(x′1，x′3)平面内，与全局坐标系中的坐标平面(x1，x3)夹角为
步骤2、确立局部坐标系下正交各向异性介质单方位反射P波走时近似公式；
(1)反射走时计算基本公式
在二维坐标系下，反射波走时T具有以下形式：
其中，x为偏移距(震源S和检波器R之间的距离)，H是水平反射面的深度(与地层的对称平面重合)。T＝T(x)表示反射P波的走时。v＝v(n)表示群速度，n为相位矢量。利用时差分析中广泛应用的表示方法：
其中，为归一化的偏移距，α为参考各向同性介质中的P波相速度，T0为参考零偏移距双程走时。公式(1)可改写为：
(2)正交各向异性介质反射走时公式
采用相速度平方近似，即v2(N)～c2(N)近似表示群速度v(N)。在二维坐标系下，正交各向异性介质的群速度可使用各向异性参数表示为：
公式(3)中，N为速度方向矢量，可以通过表示，
涉及的各向异性参数包括εx，εz，ηy，其中，εx，εz分别代表介质在坐标轴x1、x3方向P波相速度与参考速度α的差异；ηy代表介质在坐标平面(x1,x3)内P波速度的变化梯度。局部坐标系下的各向异性参数与局部坐标系下的刚度系数有关：
其中，Aij为地层介质在局部坐标系下的密度归一化刚度系数矩阵中的元素，α为参考各向同性介质中的P波相速度。
经简化，公式(3)可写为以下形式：
步骤3、基于单方位走时的局部坐标系正交各向异性介质各向异性参数反演；
(1)基于零偏移距走时计算εz
在零偏移距处，即x＝0时，公式(7)和公式(8)变为：
因此，可通过公式(12)直接求出A参数εz：
(2)构建反演行列式
根据公式(10)，公式(7)和公式(8)可改写为：
根据公式(11)构建反演行列式如下：
Gm＝d, (12)
其中，G为一N×M的系数矩阵，N为检波器的数量，M为涉及的A参数的数量，此节中M＝2。矩阵G中的行向量为以下形式：
m由2个参数组成，其形式如下：
m＝(εx，ηy)T， (14)
向量d的形式如下：
为归一化偏移距。
(3)反演行列式求解
公式(12)可以通过最小二乘的方法进行求解，反演的目标函数为：
||Gm-d||2＜κ， (16)
其中，κ为控制反演精度的参数。通过公式(17)对公式(16)进行求解：
其中，为反演算子
其中，GT为矩阵G的转置，I为N×N的单位矩阵，ε为稳定因子，ε越大，反演流程越稳定，但精度较低。
步骤4、建立全局坐标系下的各向异性参数转换公式；
在三维坐标系下，首先假设在平面上沿不同方向分布的地震测线，测线的走向角为由以测线走向为x1轴方向的局部坐标系向全局坐标系转换，转换矩阵R为：
我们将局部坐标系中的正交各向异性介质各向异性参数使用上标“′”来与全局坐标系中的正交各向异性介质各向异性参数进行区分。根据Bond变换中的6×6刚度系数矩阵转换法则，局部坐标系中的正交各向异性介质各向异性参数ε′x，η′y和ε′z可以通过以下公式转换为全局坐标系：
ε′z＝εz, (22)
局部坐标系下的各向异性参数定义见公式(6)，全局各向异性参数包括εx，εy，εz，ηx，ηy，ηz，其中，εx，εy，εz分别代表介质在坐标轴x1、x2、x3方向P波相速度与参考速度α的差异；ηx，ηy，ηz分别代表介质在坐标平面(x2,x3)、(x1,x3)、(x1,x2)内P波速度的变化梯度。全局坐标系下的各向异性参数与全局坐标系下的刚度系数有关：
其中，Aij为地层介质在全局坐标系下的密度归一化刚度系数矩阵中的元素，α为参考各向同性介质中的P波相速度。
步骤5、基于转换公式进行全局各向异性参数反演；
(1)构建反演行列式
根据公式(20)～(22)，由于沿不同方向分布的地震测线所在局部坐标系之间仅存在方位角的差异，因此各局部坐标系下的A参数εz与全局坐标系下相同，不需要进一步反演。对于剩余五个非零A参数，可通过构建以下行列式进行下一步反演。根据公式(20)和(21)，构建反演行列式如下：
G1m1＝d1, (24)
G2m2＝d2, (25)
其中，G1为N×3的系数矩阵，N为沿不同方向分布的地震测线的数量，由于公式(20)中非零A参数个数为3个，因此矩阵G1的列数为3。G2为N×2的系数矩阵，行数N为沿不同方向分布的地震测线的数量，列数为公式(21)中非零A参数的个数。矩阵G1、G2中的行向量分别为以下形式：
向量m1、m2的形式分别如下：
m1＝(εx，εy，ηz)T， (28)
m2＝(ηx，ηy)T， (29)
向量d1、d2的形式分别如下：
为了使公式(24)、(25)的反问题为适定反问题，需要系数矩阵G1和G2为正定或者超定矩阵。因此，通过多个方位地震测线所在局部坐标系下的A参数反演得到全局坐标系下的A参数，需要至少3个方位的地震反射波走时数据。
(2)反演行列式求解
公式(24)(25)可以通过最小二乘的方法进行求解，反演的目标函数为：
||G1m1-d1||2＜κ， (32)
||G2m2-d2||2＜κ， (33)
其中，κ为控制反演精度的参数。通过公式(34)对公式(32)(33)进行求解：
其中，为反演算子
其中，GT为矩阵G的转置，I为N×N的单位矩阵，ε为稳定因子，ε越大，反演流程越稳定，但精度较低。
步骤6、进行反演结果精度评价。
通过将数据协方差矩阵Cd转化为模型协方差矩阵Cm的方法评估所求各向异性参数的精度：
Cm＝HCdHT， (36)
此处模型协方差矩阵Cm表示所求的模型参数向量m的随机分布概率，若随机性不存在，则协方差矩阵为0，此时模型参数向量m没有误差。
如果准确的模型协方差矩阵通过以下方法进行计算：
Cd～σ2I， (37)
其中I为N×N的单位矩阵，σ通过χ2残差统计方法进行计算:
其中： 为公式(17)的计算结果，v为矩阵G行数与列数的差。
以下为应用本发明的几个具体实施例
实施例1
在应用本发明的一具体实施例1中，表1是本发明在数值验证过程中使用的介质模型参数，该模型包含一组6×6的刚度系数矩阵，矩阵沿对角元素对称，因此仅展示上三角矩阵。
表1实施例1中使用的各向异性模型及刚度系数矩阵表
表2是根据表1中提供的刚度系数矩阵以及参考P波速度，通过公式(23)中的各向异性参数公式计算的各向异性参数，作为精确值与反演得到的各向异性参数进行对比。同时也展示了本实施例中数据的分布范围，其中，为地震测线的走向角，x为地震测线的长度。
表2实施例1中模型的参考P波速度及各向异性参数表
图2为根据表2中数据分布范围参数绘制的实施例1中地震测线的分布范围，图3为根据公式(36)-公式(38)计算的协方差矩阵，图4为反演得到的各向异性参数与表2中精确值的对比，其误差棒根据图3中的协方差矩阵计算而来。图3、图4的结果可以证明本发明中的正交各向异性介质各向异性参数反演方法可以较精确的得到介质的各向异性参数。
实施例2
在应用本发明的一具体实施例2中，表3是本发明在数值验证过程中使用的介质模型参数，该模型包含一组6×6的刚度系数矩阵，矩阵沿对角元素对称，因此仅展示上三角矩阵。
表3实施例2中使用的各向异性模型及刚度系数矩阵表
表4是根据表3中提供的刚度系数矩阵以及参考P波速度，通过公式(23)中的各向异性参数公式计算的各向异性参数，作为精确值与反演得到的各向异性参数进行对比。同时也展示了本实施例中数据的分布范围，其中，为地震测线的走向角，x为地震测线的长度。
表4实施例2中模型的参考P波速度及各向异性参数表
图5为根据表4中数据分布范围参数绘制的实施例2中地震测线的分布范围，图6为根据公式(36)-公式(38)计算的协方差矩阵，图7为反演得到的各向异性参数与表4中精确值的对比，其误差棒根据图6中的协方差矩阵计算而来。图6、图7的结果可以证明在测线分布较为稀疏且仅有四条测线时，本发明中的正交各向异性介质各向异性参数反演方法仍然可以较精确地求得介质的各向异性参数。
实施例3
在应用本发明的一具体实施例3中，表5是本发明在数值验证过程中使用的介质模型参数，该模型包含一组6×6的刚度系数矩阵，矩阵沿对角元素对称，因此仅展示上三角矩阵。
表5实施例3中使用的各向异性模型及刚度系数矩阵表
表6是根据表5中提供的刚度系数矩阵以及参考P波速度，通过公式(23)中的各向异性参数公式计算的各向异性参数，作为精确值与反演得到的各向异性参数进行对比。同时也展示了本实施例中数据的分布范围，其中，为地震测线的走向角，x为地震测线的长度。
表6实施例3中模型的参考P波速度及各向异性参数表
图8为根据表6中数据分布范围参数绘制的实施例3中地震测线的分布范围，图9为根据公式(36)-公式(38)计算的协方差矩阵，图10为反演得到的各向异性参数与表6中精确值的对比，其误差棒根据图9中的协方差矩阵计算而来。图9、图10的结果可以证明在仅有部分方位角照明地情况下，本发明中的正交各向异性介质各向异性参数反演方法仍然可以较精确地求得介质的各向异性参数。
实施例4
在应用本发明的一具体实施例4中，表7是本发明在数值验证过程中使用的介质模型参数，该模型包含一组6×6的刚度系数矩阵，矩阵沿对角元素对称，因此仅展示上三角矩阵。
表7实施例4中使用的各向异性模型及刚度系数矩阵表
表8是根据表7中提供的刚度系数矩阵以及参考P波速度，通过公式(23)中的各向异性参数公式计算的各向异性参数，作为精确值与反演得到的各向异性参数进行对比。同时也展示了本实施例中数据的分布范围，其中，φ为地震测线的走向角，x为地震测线的长度。
表8实施例4中模型的参考P波速度及各向异性参数表
图11为根据表8中数据分布范围参数绘制的实施例3中地震测线的分布范围，图12为根据公式(36)-公式(38)计算的协方差矩阵，图13为反演得到的各向异性参数与表8中精确值的对比，其误差棒根据图12中的协方差矩阵计算而来。图12、图13的结果可以证明在地震测线最大偏移距较小地情况下，本发明中的正交各向异性介质各向异性参数反演方法仍然可以较精确地求得介质的各向异性参数。
最后应说明的是：以上所述仅为本发明的优选实施例而已，并不用于限制本发明，尽管参照前述实施例对本发明进行了详细的说明，对于本领域技术人员来说，其依然可以对前述实施例记载的技术方案进行修改，或者对其中部分技术特征进行等同替换。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
除说明书所述的技术特征外，均为本专业技术人员的已知技术。
Claims (16)
1.正交各向异性介质多方位反射波走时反演方法，其特征在于，该正交各向异性介质多方位反射波走时反演方法包括：
步骤1、建立全局坐标系及局部坐标系；
步骤2、确立局部坐标系下正交各向异性介质单方位反射P波走时近似公式；
步骤3、进行基于单方位走时的局部坐标系正交各向异性介质各向异性参数反演；
步骤4、建立全局坐标系下的各向异性参数转换公式；
步骤5、基于转换公式进行全局各向异性参数反演；
步骤6、进行反演结果精度评价。
2.根据权利要求1所述的正交各向异性介质多方位反射波走时反演方法，在步骤1，假设一笛卡尔坐标系作为全局坐标系，其x3轴是垂直方向的，x1，x2轴是水平方向且相互正交，轴面(x1，x2)、(x1，x3)分别与介质的两个对称面平行；背景介质为一均匀的地层，地层厚度为H，地层底部的水平反射面与(x1，x2)平面平行；在这种情况下，P波的入射下行和反射上行射线位于垂直于反射面的平面内，并且相对于反射面法线对称。
3.根据权利要求2所述的正交各向异性介质多方位反射波走时反演方法，其特征在于，在步骤1，针对地震波观测系统，定义局部坐标系，反射射线和两个使用的坐标系示意图：全局坐标系xi和局部坐标系x′i；角是震源检波器剖面的方位角，SO是入射射线，OR是反射射线。
4.根据权利要求3所述的正交各向异性介质多方位反射波走时反演方法，其特征在于，在步骤1，地震波从全局坐标系的原点处的震源S激发，在点O处反射，传播至检波器R；将点S和R所在的轴定义为局部坐标系的x′1轴，与x′1轴在同一水平面且与x′1轴垂直的轴为x′2轴，垂直于x′1轴、x′2轴的轴为x′3轴；局部坐标系的原点与震源S一致，并且还表示全局坐标系xi的原点，垂直轴x3在局部和全局两个坐标系中也是相同的；反射射线位于局部坐标系的(x′1，x′3)平面内，与全局坐标系中的坐标平面(x1，x3)夹角为
5.根据权利要求1所述的正交各向异性介质多方位反射波走时反演方法，其特征在于，在步骤2，在二维坐标系下，反射波走时T具有以下形式：
其中，x为偏移距即震源S和检波器R之间的距离，H是水平反射面的深度，与地层的对称平面重合；T＝T(x)表示反射P波的走时；v＝v(n)表示群速度，n为相位矢量；利用时差分析中广泛应用的表示方法：
其中，为归一化的偏移距，α为参考各向同性介质中的P波相速度，T0为参考零偏移距双程走时；公式(1)可改写为：
6.根据权利要求5所述的正交各向异性介质多方位反射波走时反演方法，其特征在于，在步骤2，采用相速度平方近似，即v2(N)～c2(N)近似表示群速度v(N)；在二维坐标系下，正交各向异性介质的群速度可使用各向异性参数表示为：
公式(3)中，N为速度方向矢量，可以通过表示，
涉及的各向异性参数包括εx，εz，ηy，其中，εx，εz分别代表介质在坐标轴x1、x3方向P波相速度与参考速度α的差异；ηy代表介质在坐标平面(x1,x3)内P波速度的变化梯度；局部坐标系下的各向异性参数与局部坐标系下的刚度系数有关：
其中，Aij为地层介质在局部坐标系下的密度归一化刚度系数矩阵中的元素，α为参考各向同性介质中的P波相速度；
经简化，公式(3)可写为以下形式：
7.根据权利要求6所述的正交各向异性介质多方位反射波走时反演方法，其特征在于，步骤3包括：
步骤31，基于零偏移距走时计算εz；
步骤32，构建反演行列式；
步骤33，求解反演行列式。
8.根据权利要求7所述的正交各向异性介质多方位反射波走时反演方法，其特征在于，在步骤31，在零偏移距处，即x＝0时，公式(7)和公式(8)变为：
因此，可通过公式(12)直接求出A参数εz：
9.根据权利要求8所述的正交各向异性介质多方位反射波走时反演方法，其特征在于，在步骤32，根据公式(10)，公式(7)和公式(8)可改写为：
根据公式(11)构建反演行列式如下：
Gm＝d, (12)
其中，G为一N×M的系数矩阵，N为检波器的数量，M为涉及的A参数的数量，此节中M＝2；矩阵G中的行向量为以下形式：
m由2个参数组成，其形式如下：
m＝(εx，ηy)T， (14)
向量d的形式如下：
为归一化偏移距，Tobs为检波器接收到的走时。
10.根据权利要求9所述的正交各向异性介质多方位反射波走时反演方法，其特征在于，在步骤33，公式(12)可以通过最小二乘的方法进行求解，反演的目标函数为：
||Gm-d||2＜κ， (16)
其中，κ为控制反演精度的参数；通过公式(17)对公式(16)进行求解：
其中，为反演算子
其中，GT为矩阵G的转置，I为N×N的单位矩阵，ε为稳定因子，ε越大，反演流程越稳定，但精度较低。
11.根据权利要求10所述的正交各向异性介质多方位反射波走时反演方法，其特征在于，在步骤4，在三维坐标系下，首先假设在平面上沿不同方向分布的地震测线，测线的走向角为由以测线走向为x1轴方向的局部坐标系向全局坐标系转换，转换矩阵R为：
将局部坐标系中的正交各向异性介质各向异性参数使用上标′来与全局坐标系中的正交各向异性介质各向异性参数进行区分；根据Bond变换中的6×6刚度系数矩阵转换法则，局部坐标系中的正交各向异性介质各向异性参数ε′x，η′y和ε′z可以通过以下公式转换为全局坐标系：
ε′z＝εz, (22)
局部坐标系下的各向异性参数定义见公式(6)，全局各向异性参数包括εx，εy，εz，ηx，ηy，ηz，其中，εx，εy，εz分别代表介质在坐标轴x1、x2、x3方向P波相速度与参考速度α的差异；ηx，ηy，ηz分别代表介质在坐标平面(x2,x3)、(x1,x3)、(x1,x2)内P波速度的变化梯度；全局坐标系下的各向异性参数与全局坐标系下的刚度系数有关：
其中，Aij为地层介质在全局坐标系下的密度归一化刚度系数矩阵中的元素，α为参考各向同性介质中的P波相速度。
12.根据权利要求11所述的正交各向异性介质多方位反射波走时反演方法，其特征在于，步骤5包括：
步骤51，构建反演行列式；
步骤52，进行反演行列式求解。
13.根据权利要求12所述的正交各向异性介质多方位反射波走时反演方法，其特征在于，在步骤51，根据公式(20)～(22)，由于沿不同方向分布的地震测线所在局部坐标系之间仅存在方位角的差异，因此各局部坐标系下的A参数εz与全局坐标系下相同，不需要进一步反演；对于剩余五个非零A参数，可通过构建以下行列式进行下一步反演；根据公式(20)和(21)，构建反演行列式如下：
G1m1＝d1, (24)
G2m2＝d2, (25)
其中，G1为N×3的系数矩阵，N为沿不同方向分布的地震测线的数量，由于公式(20)中非零A参数个数为3个，因此矩阵G1的列数为3；G2为N×2的系数矩阵，行数N为沿不同方向分布的地震测线的数量，列数为公式(21)中非零A参数的个数；矩阵G1、G2中的行向量分别为以下形式：
向量m1、m2的形式分别如下：
m1＝(εx，εy，ηz)T， (28)
m2＝(ηx，ηy)T， (29)
向量d1、d2的形式分别如下：
为了使公式(24)、(25)的反问题为适定反问题，需要系数矩阵G1和G2为正定或者超定矩阵；因此，通过多个方位地震测线所在局部坐标系下的A参数反演得到全局坐标系下的A参数，需要至少3个方位的地震反射波走时数据。
14.根据权利要求13所述的正交各向异性介质多方位反射波走时反演方法，其特征在于，在步骤52，公式(24)(25)可以通过最小二乘的方法进行求解，反演的目标函数为：
||G1m1-d1||2＜κ， (32)
||G2m2-d2||2＜κ， (33)
其中，κ为控制反演精度的参数；通过公式(34)对公式(32)(33)进行求解：
其中，为反演算子
其中，GT为矩阵G的转置，I为N×N的单位矩阵，ε为稳定因子，ε越大，反演流程越稳定，但精度较低。
15.根据权利要求14所述的正交各向异性介质多方位反射波走时反演方法，其特征在于，在步骤6，通过将数据协方差矩阵Cd转化为模型协方差矩阵Cm的方法评估所求各向异性参数的精度：
Cm＝HCdHT， (36)
此处模型协方差矩阵Cm表示所求的模型参数向量m的随机分布概率，若随机性不存在，则协方差矩阵为0，此时模型参数向量m没有误差，HT为矩阵H的转置。
16.根据权利要求15所述的正交各向异性介质多方位反射波走时反演方法，其特征在于，在步骤6，如果准确的模型协方差矩阵通过以下方法进行计算：
Cd～σ2I， (37)
其中I为N×N的单位矩阵，σ通过χ2残差统计方法进行计算:
其中： 为公式(17)的计算结果，v为矩阵G行数与列数的差。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310595476.3A CN116755148A (zh) | 2023-05-25 | 2023-05-25 | 正交各向异性介质多方位反射波走时反演方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310595476.3A CN116755148A (zh) | 2023-05-25 | 2023-05-25 | 正交各向异性介质多方位反射波走时反演方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN116755148A true CN116755148A (zh) | 2023-09-15 |
Family
ID=87948713
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310595476.3A Pending CN116755148A (zh) | 2023-05-25 | 2023-05-25 | 正交各向异性介质多方位反射波走时反演方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116755148A (zh) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US11953633B1 (en) * | 2022-10-24 | 2024-04-09 | China University Of Petroleum (East China) | Method, device and computer device for decoupling anisotropic elastic wave |
-
2023
- 2023-05-25 CN CN202310595476.3A patent/CN116755148A/zh active Pending
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US11953633B1 (en) * | 2022-10-24 | 2024-04-09 | China University Of Petroleum (East China) | Method, device and computer device for decoupling anisotropic elastic wave |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109425896B (zh) | 白云岩油气储层分布预测方法及装置 | |
| CN101329405B (zh) | 一种简单的多参数地震反演方法 | |
| Pei et al. | Velocity calibration for microseismic monitoring: A very fast simulated annealing (VFSA) approach for joint-objective optimization | |
| CN101551466B (zh) | 一种利用与偏移距有关的地震属性提高油气储层预测精度的方法 | |
| CN110133715B (zh) | 一种基于初至时差和波形叠加的微地震震源定位方法 | |
| WO2017167191A1 (zh) | 地震数据处理方法和装置 | |
| CN106094029A (zh) | 利用偏移距矢量片地震数据预测储层裂缝的方法 | |
| CN106249295B (zh) | 一种井中微地震p、s波联合快速定位方法及系统 | |
| CN106556861B (zh) | 一种基于全方位地震资料的方位avo反演方法 | |
| CN110231652B (zh) | 一种基于密度的含噪声应用空间聚类的地震相提取方法 | |
| CN102053261A (zh) | 一种地震数据处理方法 | |
| CN108957548A (zh) | 一种多波多分量联合观测地震页岩气富集区预测技术 | |
| Thiel et al. | Comparison of acoustic and elastic full‐waveform inversion of 2D towed‐streamer data in the presence of salt | |
| CN116755148A (zh) | 正交各向异性介质多方位反射波走时反演方法 | |
| US20190346581A1 (en) | Methods for determining transversely isotropic-elastic constants from borehole sonic velocities in strongly transversely-isotropic formations | |
| Wang et al. | Analysis and estimation of an inclusion-based effective fluid modulus for tight gas-bearing sandstone reservoirs | |
| CN103558637B (zh) | 基于三分量传感器的远探测方法 | |
| Pan et al. | Detection of natural tilted fractures from azimuthal seismic amplitude data based on linear-slip theory | |
| CN109696704A (zh) | 一种基于纵波阻抗约束的地震各向异性δ建模方法 | |
| Guo et al. | Becoming effective velocity-model builders and depth imagers, Part 2—The basics of velocity-model building, examples and discussions | |
| CN109324344A (zh) | 基于纯纵波和拟声波反演的页岩厚度预测方法及系统 | |
| CN113671566B (zh) | 一种基于深度域地震数据计算裂缝参数的方法 | |
| WO2017015954A1 (zh) | 一种地震信号处理方法、装置和系统 | |
| RU2705519C2 (ru) | Способ получения мигрированных сейсмических изображений геологических сред по данным сейсморазведки 2d | |
| CN107942373B (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 |