CN113671574B - 基于频域的最小二乘共轭梯度迭代的地震正演模拟方法 - Google Patents
基于频域的最小二乘共轭梯度迭代的地震正演模拟方法 Download PDFInfo
- Publication number
- CN113671574B CN113671574B CN202110950838.7A CN202110950838A CN113671574B CN 113671574 B CN113671574 B CN 113671574B CN 202110950838 A CN202110950838 A CN 202110950838A CN 113671574 B CN113671574 B CN 113671574B
- Authority
- CN
- China
- Prior art keywords
- iteration
- model
- forward modeling
- conjugate gradient
- matrix
- 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 87
- 238000004364 calculation method Methods 0.000 claims abstract description 37
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 10
- 230000009466 transformation Effects 0.000 claims abstract description 7
- 239000011159 matrix material Substances 0.000 claims description 38
- 238000004088 simulation Methods 0.000 claims description 17
- 230000008569 process Effects 0.000 claims description 12
- 238000007781 pre-processing Methods 0.000 claims description 9
- 238000010521 absorption reaction Methods 0.000 claims description 8
- 238000012952 Resampling Methods 0.000 claims description 5
- 230000009191 jumping Effects 0.000 claims description 3
- OIGNJSKKLXVSLS-VWUMJDOOSA-N prednisolone Chemical compound O=C1C=C[C@]2(C)[C@H]3[C@@H](O)C[C@](C)([C@@](CC4)(O)C(=O)CO)[C@@H]4[C@@H]3CCC2=C1 OIGNJSKKLXVSLS-VWUMJDOOSA-N 0.000 claims description 2
- 238000010586 diagram Methods 0.000 description 7
- 230000008901 benefit Effects 0.000 description 6
- 238000000354 decomposition reaction Methods 0.000 description 4
- 230000007547 defect Effects 0.000 description 2
- 239000006185 dispersion Substances 0.000 description 2
- 230000006870 function Effects 0.000 description 2
- 230000002194 synthesizing effect Effects 0.000 description 2
- BQCADISMDOOEFD-UHFFFAOYSA-N Silver Chemical compound [Ag] BQCADISMDOOEFD-UHFFFAOYSA-N 0.000 description 1
- 238000009825 accumulation Methods 0.000 description 1
- 230000006978 adaptation Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 230000004069 differentiation Effects 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000000605 extraction Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 229910052709 silver Inorganic materials 0.000 description 1
- 239000004332 silver Substances 0.000 description 1
- 230000005477 standard model 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
- 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/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
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)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开了基于频域的最小二乘共轭梯度迭代的地震正演模拟方法，所述方法包括：输入经傅里叶变换至频率域后的时谐波动方程；设置相应的正演模拟参数；根据所述正演模拟参数，通过最小二乘共轭梯度算法对所述时谐波动方程进行迭代求解。本发明的方法相较于基于LU直接求解的正演方法，可在保持同等精度的条件下减少内存消耗、提高计算效率，相较于基于BI‑CGSTAB迭代求解的正演方法，可在计算效率接近的情况下获得更宽松的使用条件和更好的计算稳定性。
Description
技术领域
本发明涉及地震勘探的技术领域。
背景技术
地震勘探作为一种基本的地球物理勘探方法，是解决油气勘探问题最有效的方法。基于波动方程的正演模拟是地震勘探的基础性工作，在地震勘探的采集、处理、解释等阶段都有着重要作用。
现有技术中的地震正演方法按照计算域可以分为时间域和频率域两类，其中时间域的算法因为实现简单，精度高得到了广泛应用，但在全波形反演的要求下，需要以正演为基础反复迭代求解波动方程，如果在时间域进行正演就会带来巨大的计算量，消耗大量的内存和时间。
相比时间域正演，频率域正演具有无时间频散、不存在稳定性问题、没有误差累积、频带选取灵活等优点，其时间域波动方程经傅里叶变换至频率域后称作时谐波动方程，该方程的求解是频率域正演的重要一环，目前常用的时谐波动方程求解方法主要有两大类：①直接求解方法，其精度高但内存消耗大，计算效率低；②迭代方法：精度相对较低，但内存消耗低，计算效率高，在三维条件下直接求解十分困难，因此迭代法的使用优势更为明显。为克服迭代法的缺陷，现有技术中提出了不同形式的迭代求解方法，其大部分都基于Krylov子空间法，如CG方法、Bi-CGSTAB方法等，但其中CG方法收敛速度较慢，优势不明显；Bi-CGSTAB方法收敛速度很快，但容易出现迭代不收敛的情况，需要较苛刻的使用条件。
如现有的专利文献CN103901472A提出的一种基于LU分解的直接求解方法，其通过建立17点格式的差分公式、构造稀疏矩阵，其后读入子波参数和速度模型，再经频率循环得到单频波场，最后对所有频率波场求反傅里叶变换得到正演结果。通过本发明实施例的方法及装置，四阶高精度正演克服了以往对dx必须等于dz的条件限制，不仅可应用于dx＝dz,也可以适用于dx≠dz的情况，且最小波长所需要的网格数在该方法中仅需要2.4个，优于其他类似方法，但该技术采用LU分解求解时谐波动方程，虽然精度高，但内存消耗大，计算效率低。
或如现有技术文献“Frequency-space domain acoustic wave simulation withthe BiCGstab,iterative method[J]”(Journal of Geophysics&Engineering,2016,13(1):70-77Du Z,Liu J,Liu W,et al.)公开了使用Bi-CGSTAB迭代方法求解大型稀疏线性方程组，其随着空间微分阶数的增加，极大地减少了内存需求和计算时间。与LU分解相比，该Bi-CGSTAB迭代的计算效率在采用四阶微分时，提高了9倍，并且其迭代的数值模拟结果等价于LU分解结果，但该方法虽然内存消耗和计算时间显著降低，但其使用条件苛刻，容易出现计算不稳定的情况。
发明内容
本发明的目在于针对现有技术的缺陷，提供一种基于频域的最小二乘共轭梯度迭代的地震正演模拟方法，其相较于基于LU直接求解的正演方法，可在保持相当精度的条件下减少内存消耗、提高计算效率，相较于基于BI-CGSTAB迭代求解的正演方法，可在计算效率接近的情况下获得更宽松的使用条件和更好的计算稳定性。
本发明的技术方案如下：
基于频域的最小二乘共轭梯度迭代的地震正演模拟方法，其包括：
1、输入经傅里叶变换至频率域后的时谐波动方程；
2、设置相应的正演模拟参数；
3、根据所述正演模拟参数，通过最小二乘共轭梯度算法对所述时谐波动方程进行迭代求解。
根据本发明的一些具体实施方式，所述模拟方法还包括：
当经过所述迭代求解获得的目标参数的计算误差收敛至预设误差的范围时，以所得目标参数进行正演模拟，输出模拟结果。
根据本发明的一些具体实施方式，所述迭代求解的过程包括：
通过差分法对所述时谐波动方程进行离散化处理，得到其线性系统模型Ap＝b(1)，其中，A表示系数矩阵，p表示压力场向量，b表示所述时谐波动方程中的震源项的离散形式，且b＝S(ω)δ(x-xs)δ(z-zs)(2)，其中，S(ω)表示傅里叶变换后频率域的震源项，δ()表示狄拉克函数，x表示横坐标自变量，z表示纵坐标自变量，xs表示震源横坐标，zs表示震源纵坐标；
通过与所述系数矩阵A具有同样稀疏结构的预处理矩阵M对所述线性系统模型进行对角预处理，得到其等价线性系统模型(AM-1)v＝b,v＝Mp(3)其中，M表示预处理矩阵，v表示中间向量；
设置待计算的目标参数x的迭代初始值x0及预设误差ε；
基于所述迭代初始值，通过所述最小二乘共轭梯度迭代算法求解所述等价线性系统模型，至其中迭代的计算误差收敛至所述预设误差ε时，获得待计算参数、即目标参数的结果，以该结果进行正演模拟。
根据本发明的一些具体实施方式，所述最小二乘共轭梯度迭代算法的求解模型如下：
其中：
sk＝αkdk(10)；
其中，k表示迭代次数，θk表示第k次迭代的最小二乘系数，表示第k次迭代时的最小二乘共轭梯度参数，/>表示第k次迭代时的第一共轭梯度参数，/>表示第k次迭代时的第二共轭梯度参数，gk＝▽f(xk)表示第k次迭代后获得的目标参数xk处的梯度矩阵，则表示xk+1处梯度矩阵gk+1的转置，yk表示前后两次迭代中的梯度差，即yk＝gk+1-gk，/>表示第k次迭代时的搜索方向矩阵dk的转置，且/>即初始搜索方向为初始值x0处的负梯度方向，其后的搜索方向通过参数/>和当前位置获得，/>表示第k+1次迭代时的最小二乘共轭梯度参数的第一搜索方向矩阵，/>表示第k+1次迭代时的最小二乘共轭梯度参数的第二搜索方向矩阵，sk表示迭代步长αk在搜索方向dk上的投影。
根据本发明的一些具体实施方式，基于所述求解模型的求解过程包括：
将目标参数的迭代初始值x0带入式(4)式获得第一次迭代的搜索方向d0和迭代步长α0；
沿着第一次迭代的搜索方向d0，使用迭代步长α0，进行对所述等价线性系统模型的第一次迭代计算，获得第一次迭代结果x1；
将所述第一次迭代结果x1带入如下的误差计算模型计算当前迭代误差ε1：
判断当前迭代误差ε1是否小于预设迭代误差ε；若小于，则输出解矩阵进行后续正演；若大于，则将所述第一次迭代结果x1作为新的迭代初始值重复其上的求解过程；
重复其上的求解过程，每次获得一个新的迭代的搜索方向dk和迭代步长αk，使用新的搜索方向和新的迭代步长计算新一轮的迭代结果xk及新的迭代误差εk，直至εk＜ε时跳出迭代，得到线性系统近似解；
以所述线性系统近似解参与正演模拟。
根据本发明的一些具体实施方式，所述差分法选自九点有限差分法。
根据本发明的一些具体实施方式，所述预设误差ε为0.01％。
根据本发明的一些具体实施方式，所述正演模拟结果包括单频数据切片、合成波场快照、地震单炮记录中的一种或多种。
根据本发明的一些具体实施方式，所述正演模拟参数包括迭代初始值x0和预设误差ε，及以下参数中的一项或多项：模型坐标系下，x方向的空间间隔dx，z方向的空间间隔dz，震源空间位置坐标，震源子波主频MF，正演模拟时长T，频率切片数量，频率间隔ΔF，PML吸收层数及吸收系数。
根据本发明的一些具体实施方式，所述速度模型选自均匀介质模型、层状介质模型、Marmousi模型、Overthrust模型中的一种或多种。
根据本发明的一些具体实施方式，所述Marmousi模型选自原始Marmousi模型或对原始Marmousi模型进行重采样后的改进Marmousi模型。
本发明具备以下有益效果：
本发明提供的基于最小二乘共轭梯度迭代的频域正演方法，能快速高效地完成频率域的正演，生成单频切片、波场快照或者地震记录，可以为波形反演提供重要帮助。
本发明的正演模拟方法对迭代条件较宽松，计算稳定，并且在保证相当精度条件下拥有较高的效率。
本发明的正演模拟方法具有可靠的模拟地震响应特征，在迭代误差0.01％的情况下与LU直接法相比，其单一频率切片计算下计算内存可减少15-20％，计算时间可减少58-88％，计算结果精确性相当，具有优异的性能。
本发明的模拟方法的模拟精度可根据实际速度模型、计算硬件环境等灵活设定，模型尺寸、频率切片数量等可根据计算设备的内存环境进行设定，适用性强、应用方便快捷。
附图说明
图1为一种具体的地震模拟正演方法的流程图；
图2为一种具体的系数矩阵结构图；
图3为一种具体的层状模型示意图；
图4为一种具体的Marmous-1模型示意图；
图5为层状模型在30Hz数据片地震正演结果示意图；
图6为层状模型在312ms合成波场快照结果示意图；
图7为Marmousi-1模型在15Hz数据片地震正演结果示意图；
图8为Marmousi-1模型在625ms合成波场快照结果示意图；
图9为两种方法下Marmousi-1模型地震单炮记录对比示意图。
具体实施方式
以下结合实施例和附图对本发明进行详细描述，但需要理解的是，所述实施例和附图仅用于对本发明进行示例性的描述，而并不能对本发明的保护范围构成任何限制。所有包含在本发明的发明宗旨范围内的合理的变换和组合均落入本发明的保护范围。
参照图1，根据本发明的技术方案，一种具体的基于最小二乘共轭梯度迭代的地震正演模拟方法包括：
S1输入经傅里叶变换至频率域后的速度模型时谐波动方程；
S2根据输入速度模型，设置相应的正演模拟参数；
S3根据设置的参数，通过最小二乘共轭梯度算法对时谐波动方程进行迭代求解；
S4判断当前迭代结果的相对误差是否满足预设误差范围，不满足则继续迭代求解；
S5在相对误差满足预设误差范围后，输出包括单频数据片、合成波场快照、或地震单炮记录的正演模拟结果。
其中，更具体的，速度模型可选用如均匀介质模型、层状介质模型等简单模型中的任一种，或Marmousi模型、Overthrust模型等标准模型中的任一种。
一些更具体的正演模拟参数可包括如：
x方向的空间间隔dx(单位：m)；
z方向的空间间隔dz(单位：m)；
震源空间位置；
震源子波主频MF(单位：Hz)；
正演模拟时长T(单位：s)；
频率切片数量；
频率间隔ΔF；
PML吸收层数及吸收系数；
待求解的未知量的迭代初始值x0及其预设误差ε。
在具体实施中，根据正演输出、如单频切片，波场快照或地震记录的需要，相应的参数设置可略作调整。
上述过程中，一种具体的最小二乘共轭梯度算法迭代求解过程包括：
S31采用差分法对所述时谐波动方程进行离散化，得到其线性系统模型如下：
Ap＝b(1)，其中A表示系数矩阵，p表示压力场向量，b表示对方程中震源项的离散，具体地，b＝S(ω)δ(x-xs)δ(z-zs)(2)，其中S(ω)表示傅里叶变换后频率域的震源项，δ()即狄拉克函数，x表示横坐标自变量，z表示纵坐标自变量，xs表示震源横坐标，zs表示震源纵坐标；
优选的，所述差分法通过现有技术中的九点有限差分法，根据其差分格式、最优化差分系数等对所述时谐波动方程进行处理，得到其离散化形式。
该优选方式所得离散后的时谐波动方程的结构为稀疏对角带状矩阵，如速度模型在x方向的网格点数为Nx，z方向的网格点数为Nz，则该矩阵是一个(Nx*Nz)*(Nx*Nz)大小的矩阵，其中非零元素最多不超过9*Nx*Nz个、0元素占绝大多数。
S32对线性系统模型Ap＝b进行预处理，得到等价的线性系统模型如下：
(AM-1)v＝b,v＝Mp (3)
其中，M表示与A具有同样稀疏结构的预处理矩阵，v表示中间向量，该等价线性系统的解p＝M-1v。
优选的，所述预处理根据所求系数矩阵具有的对角占优性质选择对角预处理，其可加速求解收敛速度。
S33设置迭代初始值x0及迭代的预设误差ε；
其中，所述初始值x0可设置为零向量，也可设置为其他值；所述预设误差ε可根据实际情况灵活选取，发明人发现，当其设置为0.01％时具有较优的综合效果，即其迭代次数、计算时间和计算精准度间较为平衡，若对计算结果的精准度要求更高，则设置的相对误差可适当降低，但同时迭代次数和计算时间会增加。
S34通过最小二乘共轭梯度迭代(Least-Squares Conjugate Gradient，LSCG)法迭代计算所述等价线性系统的解p＝M-1v；
其中，LSCG的共轭梯度参数如下：
其中，k表示迭代次数，θk表示最小二乘参数(随迭代次数k相应改变)，表示Hestenes和Stiefel在1952年提出的CG(Conjugate Gradient，共轭梯度)方法中的共轭梯度参数，如下：
其中，gk＝▽f(xk)表示xk处的梯度矩阵，则表示xk+1处梯度矩阵gk+1的转置，xk+1表示第k+1次迭代x的值，yk表示参数gk的差，yk＝gk+1-gk,/>表示搜索方向矩阵dk的转置，且即初始搜索方向为初始值x0处的负梯度方向，其后的搜索方向通过参数/>和当前位置寻找迭代步长、确定迭代方向，进而逼近线性方程的解；
表示Dai和Yuan在1999年提出的共轭梯度参数，如下：
θk由如下最小二乘问题求出：
其中，
sk＝αkdk(10)
其中，表示LSCG共轭梯度参数的搜索方向，/>表示拓展CG共轭梯度参数的搜索方向，sk表示步长αk在搜索方向dk上的投影，αk表示迭代步长。
计算求解过程可进一步包括：
S341将迭代初值x0带入式(4)式获得第一次迭代的搜索方向d0和迭代步长α0；
S342沿着第一次迭代的搜索方向d0，使用迭代步长α0进行第一次线性系统的迭代计算，获得第一次迭代的结果x1；
S343将第一次迭代的结果x1带入计算得到当前迭代误差ε1；
S344判断当前迭代误差ε1是否小于预设迭代误差ε；如果小于，则输出解矩阵进行后续正演；如果大于，则重复S341-S343，获得新的迭代的搜索方向dk和迭代步长αk，使用新的搜索方向和新的迭代步长计算新一轮的迭代结果xk，计算新的迭代误差εk，直至εk＜ε时跳出迭代，得到线性系统近似解；
S345将得到的线性系统近似解进行后续正演。
根据上述实施方式，本发明进一步提供了以下的一些仿真实施例。
以下实施例中使用到的LU直接法为通过将目标矩阵分解为一个单位下三角矩阵和一个上三角矩阵的乘积，或它们与一个置换矩阵的乘积的方式求解的方法。
实施例1
选择网格尺寸为Nx＝Nz＝100的均匀介质模型，其离散化后的系数矩阵结构如图2所示，可以看出，其系数矩阵高度稀疏，尺寸为10000*10000，从图中的局部放大部分可看到其9点差分的具体结构。
实施例2
选择如图3所示的二维层状模型，其模型网格尺寸为：Nx＝Nz＝300，空间间隔dx＝dz＝5m，设置参数包括：
正演时长1s，PML吸收层数为10，震源坐标(150，150)，震源主频30Hz，预设误差0.01％，初始向量x0为零。
根据本发明的正演模拟方法，在相对误差满足预设误差的范围后，输出的30Hz频率切片如附图5所示，312ms时刻的波场快照如附图6所示，可以看出：图5单频数据片特征清楚，上层低速层波数较下层高速层多，且二者波数关系符合速度模型理论关系速度分界面也清晰地展现出来，模拟结果正确；图6波场快照，波形特征清晰，分界面上下分别可以看到反射波透射波，模拟结果正确。
实施例3
选择如图4所示的Marmousi模型(Marmousi-1模型)，其模型网格尺寸为：Nx＝737，Nz＝751，将其横向、纵向分别都重采样1/4样点，因此模型网格尺寸变为Nx＝184，Nz＝187，空间间隔dx＝12.5m、dz＝10m，设置参数包括：
正演时长2s，PML吸收层数为50，震源坐标(50，50)，震源主频15Hz，接收点位于炮点上方z＝49处，预设误差0.01％，初始向量x0为零。
根据本发明的正演模拟方法，在相对误差满足预设误差的范围后，输出的15Hz频率切片如附图7所示，625ms时刻的波场快照如附图8所示，可以看出：图7单频数据片模型右方的高速带在数据片上有一定显示，由于提取主频较低，特征不是很清晰，模拟结果正确；图8波场快照，靠近震源的第一个强速度分界面特征清楚，投射反射波清晰，模拟结果正确。
输出的去除了PML处理区域后的地震单炮记录如附图9所示，其中，图(a)为经本发明的LSCG迭代法得到的记录，图(b)为作为对比的LU直接法得到的记录，可以看到二者的波形特征基本一致，证明本发明可得到与直接法相同的精确性。
实施例4
选择如实施例3的Marmousi-1模型，即第一模型(原始模型网格尺寸为：Nx＝737，Nz＝751，将其横向、纵向分别都重采样1/4样点，获得网格尺寸变为Nx＝184，Nz＝187的Marmousi-1模型)；
如下的Marmousi-2模型，即第二模型：
原始模型网格尺寸为：Nx＝737，Nz＝751，将其横向、纵向分别都重采样1/2样点，获得网格尺寸变为Nx＝368，Nz＝375的Marmousi-2模型；
及网格尺寸为Nx＝737、Nz＝751的原始Marmousi模型，即第三模型。
空间间隔都设置为dx＝12.5m、dz＝10m；
各模型的模拟过程在硬件参数：Intel(R)Xeon(R)Silver 4110CPU@2.1GHz(8核16线程)，RAM 32GB下单线程进行；设置参数除正演时长根据模型尺寸不同等比例延长外，其余参数与实施例3相同，同时以LU直接法作为对比。
根据本发明的正演模拟方法，在相对误差满足预设误差的范围后，各模型的计算所占内存及计算时间情况如下表所示：
Marmousi模型30Hz频率切片计算数据表
可以看出，随着模型尺寸增大，LU方法和LSCG迭代方法所占用内存均呈线性增长，但LSCG所占用内存相较LU减少了15-20％不等，而且模型越大，内存减少越多；模型尺寸增大4倍，LU方法计算时间会增加接近15倍，而LSCG迭代方法计算时间与迭代次数相关，由于其良好的收敛性，增加的计算时间并不多，相比LU方法，在第二、三模型中，计算时间减少了70％以上，第三模型甚至减少了87.98％。
以上实施例仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例。凡属于本发明思路下的技术方案均属于本发明的保护范围。应该指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下的改进和润饰，这些改进和润饰也应视为本发明的保护范围。
Claims (5)
1.基于频域的最小二乘共轭梯度迭代的地震正演模拟方法，其特征在于，包括以下步骤：
步骤S1：输入经傅里叶变换至频率域后的速度模型时谐波动方程；
步骤S2：设置相应的正演模拟参数，所述正演模拟参数包括迭代初始值x0和预设误差ε，及以下参数中的一项或多项：模型坐标系下，x方向的空间间隔dx，z方向的空间间隔dz，震源空间位置坐标，震源子波主频MF，正演模拟时长T，频率切片数量，频率间隔ΔF，PML吸收层数及吸收系数；
步骤S3：根据所述正演模拟参数，通过最小二乘共轭梯度算法对所述时谐波动方程进行迭代求解；
所述迭代求解的过程包括：
通过差分法对所述时谐波动方程进行离散化处理，得到其线性系统模型Ap＝b(1)，其中A表示系数矩阵，p表示压力场向量，b表示所述时谐波动方程中的震源项的离散形式，且b＝S(ω)δ(x-xs)δ(z-zs) (2)，其中S(ω)表示傅里叶变换后频率域的震源项，δ()表示狄拉克函数，x表示横坐标自变量，z表示纵坐标自变量，xs表示震源横坐标，zs表示震源纵坐标；
通过与所述系数矩阵A具有同样稀疏结构的预处理矩阵M对所述线性系统模型进行对角预处理，得到其等价线性系统模型(AM-1)v＝b,v＝Mp(3)，其中M表示预处理矩阵，v表示中间向量；
设置迭代初始值x0及预设误差ε；
基于所述迭代初始值，通过所述最小二乘共轭梯度算法求解所述等价线性系统模型，至其中迭代的计算误差收敛至所述预设误差ε时，获得待计算参数、即目标参数的结果，以该结果进行正演模拟；
其中，所述最小二乘共轭梯度算法的求解模型如下：
其中：
sk＝αkdk(10)；
其中k表示迭代次数，θk表示第k次迭代的最小二乘系数，表示第k次迭代时的最小二乘共轭梯度参数，/>表示第k次迭代时的第一共轭梯度参数，/>表示第k次迭代时的第二共轭梯度参数，gk＝▽f(xk)表示第k次迭代后获得的目标参数xk处的梯度矩阵，则/>表示xk+1处梯度矩阵gk+1的转置，yk表示前后两次迭代中的梯度差，即yk＝gk+1-gk，/>表示第k次迭代时的搜索方向矩阵dk的转置，且d0＝-g0,/>即初始搜索方向为初始值x0处的负梯度方向，其后的搜索方向通过参数/>和当前位置获得，/>表示第k+1次迭代时的最小二乘共轭梯度参数的第一搜索方向矩阵，/>表示第k+1次迭代时的最小二乘共轭梯度参数的第二搜索方向矩阵，sk表示迭代步长αk在搜索方向dk上的投影；
基于所述求解模型的求解过程包括：
将目标参数的迭代初始值x0带入式(4)获得第一次迭代的搜索方向d0和迭代步长α0；
沿着第一次迭代的搜索方向d0，使用迭代步长α0，进行对所述等价线性系统模型的第一次迭代计算，获得第一次迭代结果x1；
将所述第一次迭代结果x1带入如下的误差计算模型计算当前迭代误差ε1：
判断当前迭代误差ε1是否小于预设迭代误差ε；若小于，则输出解矩阵进行后续正演；若大于，则将所述第一次迭代结果x1作为新的迭代初始值重复其上的求解过程；
重复其上的求解过程，每次获得一个新的迭代的搜索方向dk和迭代步长αk，使用新的搜索方向和新的迭代步长计算新一轮的迭代结果xk及新的迭代误差εk，直至εk＜ε时跳出迭代，得到线性系统近似解；
以所述线性系统近似解参与正演模拟。
2.根据权利要求1所述的模拟方法，其特征在于：所述差分法为九点有限差分法。
3.根据权利要求1所述的模拟方法，其特征在于：所述预设误差ε为0.01％。
4.根据权利要求1所述的模拟方法，其特征在于：所述速度模型选自均匀介质模型、层状介质模型、Marmousi模型、Overthrust模型中的一种或多种。
5.根据权利要求4所述的模拟方法，其特征在于：所述Marmousi模型选自标准Marmousi-1模型或对标准Marmousi模型-1进行重采样后的改进Marmousi模型。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110950838.7A CN113671574B (zh) | 2021-08-18 | 2021-08-18 | 基于频域的最小二乘共轭梯度迭代的地震正演模拟方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110950838.7A CN113671574B (zh) | 2021-08-18 | 2021-08-18 | 基于频域的最小二乘共轭梯度迭代的地震正演模拟方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113671574A CN113671574A (zh) | 2021-11-19 |
| CN113671574B true CN113671574B (zh) | 2024-04-26 |
Family
ID=78543680
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110950838.7A Active CN113671574B (zh) | 2021-08-18 | 2021-08-18 | 基于频域的最小二乘共轭梯度迭代的地震正演模拟方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113671574B (zh) |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US5557710A (en) * | 1993-02-08 | 1996-09-17 | International Business Machines Corporation | Computer aided design system |
| CN103278848A (zh) * | 2013-04-22 | 2013-09-04 | 中山大学 | 基于mpi并行预条件迭代的地震成像正演方法 |
| CN110794457A (zh) * | 2018-08-01 | 2020-02-14 | 中国石油化工股份有限公司 | 一种层析反演的方法及系统 |
Family Cites Families (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US8274858B2 (en) * | 2009-11-12 | 2012-09-25 | Pgs Geophysical As | Method for full-bandwidth deghosting of marine seismic streamer data |
| GB2501829B (en) * | 2011-02-24 | 2019-07-17 | Chevron Usa Inc | System and method for performing reservoir simulation using preconditioning |
| US9075159B2 (en) * | 2011-06-08 | 2015-07-07 | Chevron U.S.A., Inc. | System and method for seismic data inversion |
| US10838093B2 (en) * | 2015-07-02 | 2020-11-17 | Exxonmobil Upstream Research Company | Krylov-space-based quasi-newton preconditioner for full-wavefield inversion |
-
2021
- 2021-08-18 CN CN202110950838.7A patent/CN113671574B/zh active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US5557710A (en) * | 1993-02-08 | 1996-09-17 | International Business Machines Corporation | Computer aided design system |
| CN103278848A (zh) * | 2013-04-22 | 2013-09-04 | 中山大学 | 基于mpi并行预条件迭代的地震成像正演方法 |
| CN110794457A (zh) * | 2018-08-01 | 2020-02-14 | 中国石油化工股份有限公司 | 一种层析反演的方法及系统 |
Non-Patent Citations (7)
| Title |
|---|
| 一种修正的MPS预条件共轭梯度算法;杨娇、庄杰鹏;湖南理工学院学报(自然科学版);摘要、引言、第1章 * |
| 三维VTI介质qP波正演方法研究;戚艳平;中国优秀硕士学位论文全文数据库;20090615;全文 * |
| 基于 Krylov 子空间法的频域地震数值模拟;周觅路、刘文革;第四届油气地球物理学术年会论文集;摘要、引言、第1-3章 * |
| 求解Helmholtz方程的新型稀疏近似逆预条件算法;李月卉;詹红霞;;半导体光电;20121015(05);全文 * |
| 病态总体最小二乘解算方法及应用研究;于冬冬;中国优秀硕士学位论文全文数据库;第4章 * |
| 频率域波动方程正演基于多重网格预条件的迭代算法研究;曹健;陈景波;曹书红;;地球物理学报;20150315(03);全文 * |
| 频率域波动方程高精度有限差分格式及并行模拟算法研究;李扬;中国博士学位论文全文数据库;20180115;全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113671574A (zh) | 2021-11-19 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107783190B (zh) | 一种最小二乘逆时偏移梯度更新方法 | |
| CN104977607B (zh) | 利用变步长网格声波波场模拟的时间域全波形反演方法 | |
| Liu et al. | A new kind of optimal second-order symplectic scheme for seismic wave simulations | |
| CN111766628A (zh) | 一种预条件的时间域弹性介质多参数全波形反演方法 | |
| CN112327358B (zh) | 一种粘滞性介质中声波地震数据正演模拟方法 | |
| Zhang et al. | 3D simultaneous seismic data reconstruction and noise suppression based on the curvelet transform | |
| CN108828659B (zh) | 基于傅里叶有限差分低秩分解的地震波场延拓方法及装置 | |
| CN113671574B (zh) | 基于频域的最小二乘共轭梯度迭代的地震正演模拟方法 | |
| CN111665556A (zh) | 地层声波传播速度模型构建方法 | |
| Chen et al. | A time-space domain stereo finite difference method for 3D scalar wave propagation | |
| CN109490954B (zh) | 波场正演模拟方法及装置 | |
| CN110658558A (zh) | 吸收衰减介质叠前深度逆时偏移成像方法及系统 | |
| CN114357831B (zh) | 无网格广义有限差分正演方法、装置、存储介质及设备 | |
| Shin et al. | Laplace-domain full waveform inversion using irregular finite elements for complex foothill environments | |
| Bartana* et al. | GPU implementation of minimal dispersion recursive operators for reverse time migration | |
| Song et al. | Implicit finite difference wave equation forward modeling based on partial fraction expansion | |
| Zheng et al. | Vibrator data denoising based on fractional wavelet transform | |
| Gou et al. | Complex seismic wavefield interpolation based on the Bregman iteration method in the sparse transform domain | |
| Bakulin et al. | Benchmarking 3D time-and frequency-domain solvers for FWI applications for different cluster sizes and variable number of sources | |
| Mastryukov | Optimal finite difference schemes for a wave equation | |
| CN113917526B (zh) | 基于非分裂完全匹配层吸收边界的正演模拟方法 | |
| CN115373020B (zh) | 一种基于离散小波矩量法的地震散射波场数值模拟方法 | |
| Tan et al. | 3D trapezoid-grid finite-difference time-domain (FDTD) method for seismic wave modeling | |
| CN109283589B (zh) | 一种重力场水平分量的获取方法 | |
| Liu et al. | Image fusion for travel time tomography inversion |
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 |