基于随机梯度下降回归算法的全波形反演方法与流程

1.本发明涉及地震勘探速度建模技术领域，涉及一种基于随机梯度下降法的全波形反演方法。


背景技术：

2.全波形反演在理论上比现在使用的传统走时成像技术所利用的信息更多，能达到更高的分辨率，能够通过使用全波场的信息以规避传统走时层析成像的高频近似问题。对于全波形反演技术的优化过程的实现，现有的梯度算法与似牛顿算法，结构较为简单容易实现，且收敛速度较快。但是由于优化问题存在局部最优解问题，传统算法容易落入局部极小点附近。全局梯度算法则受限于计算资源的限制，无法大规模应用。因此，研究一种能够提高全波形反演优化问题计算效率并能够避免局部极值问题的加速算法具有重要的理论价值以及实际意义。
3.在申请号：cn201410416052.7的中国专利申请中，涉及到基于梯度处理的全波形反演方法和系统。该方法包括:残差计算步骤,由初始速度模型正演得到的计算数据和野外观测数据确定数据残差；梯度计算步骤,根据数据残差在初始速度模型中的反传波场以及正演得到的正传波场求取梯度场；扰动模型计算步骤,对梯度场进行滤波得到保留地层边界信息的速度扰动模型；更新步骤,根据速度扰动模型更新初始速度模型。该发明对全波形反演中的梯度进行低通和自适应两次滤波,通过低通滤波压制随机噪声,并通过自适应滤波恢复在低通滤波过程中损失的边界信息,达到边界保真的效果。
4.在申请号：cn201310133893.2的中国专利申请中，涉及到一种频率域全波形反演地震速度建模方法，包括以下步骤：1)获取原始地震炮集记录、震源子波信息和反演所用的初始模型；2)分析步骤1)获得的信息，基于正演算法和优化算法，确定基本的反演参数和从低频到高频的全波形反演框架；3)对不同的频率，计算获得最合适的频域正反演模型网格；4)在低频反演时采用主成分分析方法，压缩参与反演的数据维度；5)判断对应不同频率的投影矩阵维度是否满足阈值转换标准，当满足该转换标准时，进入下一步，否则回到步骤4)；6)引入震源编码方法，并利用随机相位编码方法压制串扰噪声；7)判定是否满足迭代截止条件，如果满足迭代截止条件，则进入下一步；否则回到步骤6)；8)如果没有完成所有频率的反演，则回到步骤3)，直至完成所有频率的反演，获得最后的速度模型，输出速度模型。
5.在申请号：cn201710487837.7的中国专利申请中，涉及到一种基于地震记录积分的全波形反演方法，属于地震勘探的反演地下介质速度模型技术领域。本发明包括如下步骤：获取观测地震记录、去噪处理，提取地震子波；计算模拟地震记录；对去噪后的观测地震记录和模拟地震记录分别进行时间域积分；将积分后的两组数据作差，并将之作为反传震源，计算反传波场；计算正传波场；计算梯度；计算共轭梯度；根据线性搜索方法，计算迭代更新的步长；计算模型的修正量，并更新速度模型；判断更新后的模型是否满足给定的条件，如果满足则输出结果，否则，将更新后的速度模型作为新的初始速度模型，重新开始计算。解决了全波形反演由于缺少低频地震数据而使反演陷入局部极值，导致反演失败的问
题。
6.以上现有技术均与本发明有较大区别，未能解决我们想要解决的技术问题，为此我们发明了一种新的基于随机梯度下降回归算法的全波形反演方法。


技术实现要素：

7.本发明的目的是提供一种可以提高全波形反演方法精度，解决局部极值问题的基于随机梯度下降法的全波形反演方法。
8.本发明的目的可通过如下技术措施来实现：基于随机梯度下降回归算法的全波形反演方法，该基于随机梯度下降回归算法的全波形反演方法包括：
9.步骤1，确定初始模型；
10.步骤2，对模型进行正演分析；
11.步骤3，使用正演后得到的模拟数据和观测数据的差异计算目标函数；
12.步骤4，对正演结果和观测结果之间定义的目标函数使用伴随状态法构造梯度；
13.步骤5，借由得到的梯度使用随机梯度下降法对初始模型进行优化；
14.步骤6，重复步骤2、3、4、5直到模型收敛。
15.本发明的目的还可通过如下技术措施来实现：
16.在步骤1中，取得波包反演之后得到的模型以作为初始模型和原始地震数据。
17.在步骤1中，对数据进行预处理，包括筛选信道，数据变换，噪声的压制。
18.在步骤2中，对初始模型正演，得到正演地震记录或声波速度波场剖面；在正演模拟中，首先对声波方程进行网格化，采用有限差分方法将模型划分为规则网格，然后利用震源信息和边界条件得到正演结果的地震波场，并且就实际观测波形采集观测波场，对观测数据进行噪声处理以及数据变换预处理。
19.在步骤2中，对模型进行正演，需要一个给定的震源和波动方程，对于声波方程，定义一个声波方程：
[0020][0021]
其中t为时间，dt为时间的微分，x、y为场的任一点的横纵坐标，x0,y0为场的给定一点的横纵坐标，u为位移，du为位移的微分，c为波速，δ为拉普拉斯算子，q为一个位于(x0,y0)处的震源，同时给定区域边界，使用有限差分法对计算波场u可以得到正演结果d
pre
。
[0022]
在步骤3中，利用步骤2中计算得到的正演波场d
pre
，让其和观测数据d
obs
作差，然后计算数据残差的l2范数，即
[0023][0024]
其中ne为每个单元e。
[0025]
在步骤4中，基于伴随状态法原理，将数据残差作为虚拟震源在检波器处进行反传模拟以得到新的波场，将这一反传得到的反传波场和正传波场进行互相关运算，由于伴随波场的特性，将其互相关叠加后可以得到目标函数的梯度g。
[0026]
在步骤4中，伴随状态法构造的伴随源产生的波场为v，则梯度g可由下式求出：
[0027][0028]
其中m为反演的模型，f(mk)表示反演目标函数，t为时间，u为位移,x表示空间坐标，mk指迭代第k步的m，t为整个周期，s表示整个空间上,v为伴随状态法构造的伴随源产生的波场。
[0029]
在步骤5中，使用随机梯度下降法，利用按照步骤4中得到的k次迭代后的残差信息以及目标函数的梯度g来计算下一次模型的更新方向，并利用线性搜索方法确定其迭代的下降步长。
[0030]
在步骤5中，基于随机梯度下降法的算法原理更新模型迭代方向的流程为：
[0031]
对于目标函数f(m)，在某点ms处按二阶泰勒展开，
[0032][0033]
其中h为海森矩阵,g为目标函数的梯度，g
t
为梯度矩阵的转置，梯度下降
[0034]ms+1
＝m
s-lk*g
ꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)
[0035]
其中ms指对m做的第s次迭代后的m，将(5)代入(4)计算可得目标函数的变化为
[0036][0037]
对于随机梯度下降法，随机选取一个样本进行目标函数的更新则有
[0038]mk+1
＝m
k-lk*gkꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(7)
[0039]
其中，k为迭代次数，lk和gk分别表示第k次迭代的加权系数和目标函数的梯度。
[0040]
在步骤5中，对于加权系数lk使用线性搜索方法确定，这一系数也就是在模型迭代方向上的最佳步长，采用二分法搜索这些线性搜索方法计算得到该系数。
[0041]
在步骤6中，在设定完终止条件后，对于步骤5所得到的更新方向更新模型后进行线判断接受或放弃新解，最后判断终止条件，满足条件则停止迭代，否则更新模型重新进行步骤2至步骤5。
[0042]
本发明针对地震全波形反演容易落入局部极小点附近的问题，提供了一种精确的随机梯度下降算法反演方法。在常规的全波形反演算法中引入随机梯度下降法替代拟牛顿法，即在每一次迭代中随机选取样本进行梯度下降。与现有技术相比，本发明的有益效果是：
[0043]
本发明利用了两次迭代更新的参数模型进行了一种随机取样的梯度计算方法，计算得到了最接近全局最优解的参数模型，使下次迭代的方向与步长得到确定。在每次计算中，相比传统梯度下降方法随机梯度下降大幅减少了计算量，提高了计算效率。相比牛顿类方法，本发明采用的随机梯度下降法避免了牛顿类方法容易落入局部最小值的问题，更容易达到全局最优解。在相同精度下，比传统梯度下降有效地提高了全波形反演流程的计算效率。
附图说明
[0044]
图1为本发明的基于随机梯度下降回归算法的全波形反演方法的一具体实施例的流程图；
[0045]
图2为本发明的一具体实施例中基于随机梯度下降回归算法的全波形反演初始模型的示意图；
[0046]
图3为本发明的一具体实施例中基于随机梯度下降回归算法的全波形反演方法迭代中结果的示意图；
[0047]
图4为本发明的一具体实施例中基于随机梯度下降回归算法的全波形反演方法迭代完成的结果示意图；
[0048]
图5为本发明的一实施例中基于随机梯度下降回归算法的全波形反演方法和常规数学梯度下降相关方法的效果对比图。
具体实施方式
[0049]
为使本发明的上述和其他目的、特征和优点能更明显易懂，下文特举出较佳实施例，并配合附图所示，作详细说明如下。具体实施例和应用中所提供的描述信息仅用于示例解释本发明，并不能作为本发明的限定。本发明所定义的一般原则，在不脱离本发明的范围下可以在其他实施例中适用。
[0050]
本发明的基于随机梯度下降回归算法的全波形反演方法，首先确定初始模型，取得基于地震实际观测数据的观测波场；然后对初始模型进行正演分析，得到正传波场；然后测算正演后得到的模拟数据和观测数据的差值，计算目标函数；之后，对正演结果和观测结果之间定义的目标函数使用伴随状态法构造梯度；步骤5：借由得到的梯度使用随机梯度下降法计算参数模型更新的迭代方向，并采用线性搜索的方法更新步长，对模型进行更新；最后判断模型是否符合收敛条件，如果不符合，则继续重复上述步骤直到模型收敛。
[0051]
以下为应用本发明的几个具体实施例
[0052]
实施例1
[0053]
在应用本发明的一具体实施例1中，该基于随机梯度下降法的全波形反演目标函数优化算法包括：
[0054]
步骤1，确定初始模型；取得为了进行全波形反演问题要使用的初始地质模型，通常为层状模型。
[0055]
步骤2，对模型进行正演分析；对初始模型正演，得到正演地震记录或声波速度波场剖面；在正演模拟中，首先对声波方程进行网格化，采用有限差分方法将模型划分为规则网格，然后利用震源信息和边界条件得到正演结果的地震波场。并且就实际观测波形采集观测波场。并对观测数据进行噪声处理以及数据变换预处理。
[0056]
步骤3，使用正演后得到的模拟数据和观测数据的差异计算目标函数；
[0057]
利用步骤2中计算得到的正演波场d
pre
，让其和观测数据d
obs
作差，然后计算数据残差的l2范数，即
[0058][0059]
步骤4，对正演结果和观测结果之间定义的目标函数使用伴随状态法构造梯度；基于伴随状态法原理，将数据残差作为虚拟震源在检波器处进行反传模拟以得到新的波场，将这一反传得到的反传波场和正传波场进行互相关运算，由于伴随波场的特性，将其互相关叠加后可以得到目标函数的梯度。
[0060]
步骤5，借由得到的梯度使用随机梯度下降法对初始模型进行优化；使用随机梯度下降法，利用按照步骤4中得到的k次迭代后的残差信息以及目标函数的梯度来计算下一次模型的更新方向。并利用线性搜索等方法确定其迭代的下降步长。
[0061]
步骤6，重复步骤2、3、4、5直到模型收敛。在设定完终止条件后，对于步骤5所得到的更新方向更新模型后进行线判断接受或放弃新解，最后判断终止条件，满足条件则停止迭代，否则更新模型重新进行步骤2至步骤5。
[0062]
实施例2
[0063]
在应用本发明的一具体实施例2中，如图1所示，图1为本发明的基于随机梯度下降回归算法的全波形反演方法的流程图；该基于随机梯度下降回归算法的全波形反演方法包括：
[0064]
在步骤101，取得波包反演之后得到的模型以作为初始模型和原始地震数据。对数据进行预处理，其中包括筛选信道，数据变换，噪声的压制等方法。
[0065]
在步骤102中，对初始模型进行正演分析，为了求取整个模型的波场，将模型进行规则网格划分，之后使用有限差分法，使用数值模拟方法对模型进行正演模拟，所得到的波场被用作正演波场d
pre
。
[0066]
在步骤102中，对模型进行正演，需要一个给定的震源和波动方程，以声波方程为例，定义一个声波方程
[0067][0068]
其中c为波速，δ为拉普拉斯算子，q为一个位于(x0,y0)处的震源，同时给定区域边界。使用有限差分法对计算波场u可以得到正演结果d
pre
。
[0069]
在步骤103中，利用步骤2中计算得到的正演波场d
pre
，让其和观测数据d
obs
作差，然后计算数据残差的l2范数，即
[0070][0071]
在步骤104中，首先依据伴随状态法原理，将数据残差作为伴随震源将其放置在检波器上，在模型中进行共轭反传，得到模拟的反传波场，然后再将在步骤2中得到的正传波场和反传波场进行互相关运算，利用伴随波场的特性可以通过互相关叠加的方式得到目标函数的梯度g。
[0072]
其中伴随状态法构造的伴随源产生的波场为v，则梯度g可由下式求出：
[0073][0074]
在步骤105中，使用随机梯度下降法，利用按照步骤4中得到的k次迭代后的残差信息以及目标函数的梯度g来计算下一次模型的更新方向。并利用线性搜索等方法确定其迭代的下降步长。
[0075]
在步骤105中，基于随机梯度下降法的算法原理更新模型迭代方向的流程为：
[0076]
对于目标函数f(m)，在某点ms处按二阶泰勒展开，
[0077][0078]
其中h为海森矩阵，梯度下降
[0079]ms+1
＝m
s-lk*g
ꢀꢀꢀ
(5)
[0080]
将(5)代入(4)计算可得目标函数的变化为
[0081][0082]
对于随机梯度下降法，随机选取一个样本进行目标函数的更新则有
[0083]mk+1
＝m
k-lk*gkꢀꢀꢀꢀ
(7)
[0084]
其中，k为迭代次数，lk和gk分别表示第k次迭代的加权系数和目标函数的梯度。
[0085]
在步骤105中，对于加权系数lk使用线性搜索方法确定，这一系数也就是在模型迭代方向上的最佳步长，采用二分法搜索等线性搜索方法计算得到该系数。
[0086]
在步骤106中，在设定完终止条件后，对于步骤5所得到的更新方向更新模型后进行线判断接受或放弃新解，最后判断终止条件，满足条件则停止迭代，否则更新模型重新进行步骤102至步骤105。当停止迭代后，输出更新后的模型，流程结束。
[0087]
实施例3
[0088]
在应用本发明的一具体实施例3中，如图2-4所示，图2-4为采用本发明的一实施例，在使用基于随机梯度下降法的全波形反演目标函数优化算法对全波形反演的结果。可以看出，在全波形反演的过程中，随机梯度下降法使得模型从初始模型更为精确地迭代至全局最优解，而且剖面的构造清晰，有利于下一步分析。与传统牛顿类方法使用逆海森矩阵快速下降不同，随机梯度下降法由于取样的随机性，不易像牛顿法一般落入局部极小值中，使迭代无法继续进行。随机梯度下降法利用了多个样本特征，能够绕过下降过程中的鞍点，精确表达图像的边缘以及特征
[0089]
图5为本发明的一实施例中对于普通梯度下降法和随机梯度下降法的优化算法的效果对比图。在处理数据的情况下，由于随机梯度下降法不用在每一步中都对整个网格的梯度进行计算，而是只计算一个样本，所以大大减少了计算量，所求解的矩阵也比传统方法更简单。同时由于取样是随机性，也保有了梯度下降法所具有的准确性。
[0090]
最后应说明的是：以上所述仅为本发明的优选实施例而已，并不用于限制本发明，尽管参照前述实施例对本发明进行了详细的说明，对于本领域技术人员来说，其依然可以对前述实施例记载的技术方案进行修改，或者对其中部分技术特征进行等同替换。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
[0091]
除说明书所述的技术特征外，均为本专业技术人员的已知技术。

技术特征：
1.基于随机梯度下降法的全波形反演目标函数优化算法，其特征在于，该基于随机梯度下降法的全波形反演目标函数优化算法包括：步骤1，确定初始模型；步骤2，对模型进行正演分析；步骤3，使用正演后得到的模拟数据和观测数据的差异计算目标函数；步骤4，对正演结果和观测结果之间定义的目标函数使用伴随状态法构造梯度；步骤5，借由得到的梯度使用随机梯度下降法对初始模型进行优化；步骤6，重复步骤2、3、4、5直到模型收敛。2.根据权利要求1所述的基于随机梯度下降法的全波形反演目标函数优化方法，其特征在于，在步骤1中，取得波包反演之后得到的模型以作为初始模型和原始地震数据。3.根据权利要求2所述的基于随机梯度下降法的全波形反演目标函数优化方法，其特征在于，在步骤1中，对数据进行预处理，包括筛选信道，数据变换，噪声的压制。4.根据权利要求1所述的基于随机梯度下降法的全波形反演目标函数优化算法，其特征在于，在步骤2中，对初始模型正演，得到正演地震记录或声波速度波场剖面；在正演模拟中，首先对声波方程进行网格化，采用有限差分方法将模型划分为规则网格，然后利用震源信息和边界条件得到正演结果的地震波场，并且就实际观测波形采集观测波场，对观测数据进行噪声处理以及数据变换预处理。5.根据权利要求4所述的基于随机梯度下降法的全波形反演目标函数优化算法，其特征在于，在步骤2中，对模型进行正演，需要一个给定的震源和波动方程，对于声波方程，定义一个声波方程：其中t为时间，dt为时间的微分，x、y为场的任一点的横纵坐标，x0,y0为场的给定一点的横纵坐标，u为位移，du为位移的微分，c为波速，δ为拉普拉斯算子，q为一个位于(x0,y0)处的震源，同时给定区域边界，使用有限差分法对计算波场u可以得到正演结果d
pre
。6.根据权利要求5所述的基于随机梯度下降法的全波形反演目标函数优化算法的目标函数计算，其特征在于，在步骤3中，将正演得到的速度波场数据和观测数据作差以得到数据残差，并且将所得到的残差的l2范数作为反演的目标函数，即其中，f(m)表示反演目标函数，d
pre
表示正演模拟数据，d
obs
表示观测数据。7.根据权利要求1所述的基于随机梯度下降法的全波形反演目标函数优化算法的梯度计算，其特征在于，在步骤4中，基于伴随状态法原理，将数据残差作为虚拟震源在检波器处进行反传模拟以得到新的波场，将这一反传得到的反传波场和正传波场进行互相关运算，由于伴随波场的特性，将其互相关叠加后可以得到目标函数的梯度g。8.根据权利要求7所述的基于随机梯度下降法的全波形反演目标函数优化算法的梯度计算，其特征在于，在步骤4中，伴随状态法构造的伴随源产生的波场为v，则梯度g可由下式
求出：其中m为反演的模型，f(m
k
)表示反演目标函数，t为时间，u为位移,x表示空间坐标，m
k
指迭代第k步的m，t为周期，s表示整个空间上,v为伴随状态法构造的伴随源产生的波场。9.根据权利要求1所述的基于随机梯度下降法的全波形反演目标函数优化算法，其特征在于，在步骤5中，使用随机梯度下降法，利用按照步骤4中得到的k次迭代后的残差信息以及目标函数的梯度g来计算下一次模型的更新方向，并利用线性搜索方法确定其迭代的下降步长。10.根据权利要求9所述的基于随机梯度下降法的全波形反演目标函数优化算法，其特征在于，在步骤5中，基于随机梯度下降法的算法原理更新模型迭代方向的流程为：对于目标函数f(m)，在某点m
s
处按二阶泰勒展开，其中f(m)表示反演目标函数，h为海森矩阵，g为目标函数的梯度，g
t
为梯度矩阵的转置,梯度下降m
s+1
＝m
s-l
k
*g
ꢀꢀꢀꢀ
(5)其中m
s
指对m做的第s次迭代后的m，将(5)代入(4)计算可得目标函数的变化为对于随机梯度下降法，随机选取一个样本进行目标函数的更新则有m
k+1
＝m
k-l
k
*g
k
ꢀꢀꢀꢀ
(7)其中，k为迭代次数，l
k
和g
k
分别表示第k次迭代的加权系数和目标函数的梯度。11.根据权利要求10所述的基于随机梯度下降回归算法的全波形反演方法，其特征在于，在步骤5中，对于加权系数l
k
使用线性搜索方法确定，这一系数也就是在模型迭代方向上的最佳步长，采用二分法搜索这些线性搜索方法计算得到该系数。12.根据权利要求1所述的基于随机梯度下降回归算法的全波形反演方法，其特征在于，在步骤6中，在设定完终止条件后，对于步骤5所得到的更新方向更新模型后进行线判断接受或放弃新解，最后判断终止条件，满足条件则停止迭代，否则更新模型重新进行步骤2至步骤5。

技术总结
本发明提供一种基于随机梯度下降法的全波形反演目标函数优化算法，包括：步骤1,确定初始模型；步骤2，对初始模型进行正演分析；步骤3：测算正演后得到的模拟数据和观测数据的差值，计算目标函数；步骤4：对正演结果和观测结果之间定义的目标函数使用伴随状态法构造梯度；步骤5：借由得到的梯度使用随机梯度下降法对初始模型进行优化；步骤6：重复步骤2、3、4、5直到模型收敛。该基于随机梯度下降法的全波形反演优化算法借助了随机梯度下降法通过使用随机梯度下降优化算法对全局最优解的搜索能力，提高了全波形反演优化过程的精度，提高了对局部极值问题的抵抗能力。了对局部极值问题的抵抗能力。了对局部极值问题的抵抗能力。


技术研发人员：曲志鹏 朱剑兵 王兴谋 梁鸿贤 姜蕾 徐宏斌 刘鸽
受保护的技术使用者：中国石油化工股份有限公司胜利油田分公司物探研究院
技术研发日：2022.03.23
技术公布日：2023/10/6