CN118033743A - 一种基于自然分频时频域能量谱的弹性波反演方法 - Google Patents

Abstract

本发明涉及一种基于自然分频时频域能量谱的弹性波反演方法，本方法利用S变换提取地震数据的单频能量谱信息并嵌套频率点循环实现自然分频多尺度反演，在低频缺失情况下建立更为精确的纵横波速度初始模型，并在此基础上利用多尺度策略进行弹性波全波形反演的分频组迭代更新，获得精确的纵横波速度模型，提高低频缺失条件下弹性波全波形反演的结果精确度和可靠性。

Description

一种基于自然分频时频域能量谱的弹性波反演方法

技术领域

本发明涉及地震勘探反演技术领域，具体地，涉及一种基于自然分频时频域能量谱的弹性波反演方法。

背景技术

速度建模问题是地震勘探的最突出问题之一。速度建模的精度会影响地震成像的真实性和可靠性，干扰地震资料的解释。如何建立高精度的速度模型是解决问题的关键。速度建模方法主要有三类：一类是走时层析，其仅利用地震记录的走时信息，只能得到精度较低的速度模型；第二类是速度分析，这种建模方法受操作人员主观因素影响；第三类是全波形反演方法，这种方法能够利用了地震记录中的全部信息，可以获得高精度的速度模型。

弹性波全波形反演相比于声波全波形反演非线性更强。弹性波全波形反演的难点在于：横波速度模型反演对初始模型的依赖性更强和进行弹性波反演时的巨大计算量，前者影响到弹性波全波形反演中横波速度模型的准确性，而后者则直接关系到相关技术能够应用于生产。

Sirgue提出了一种基于时间域正演，频率域反演的混合域全波形反演方法，该方法结合了时间域正演和频率域反演的优势，通过在时间域正演过程中加入傅里叶变换以获取频率域波场，相比于频率域全波形反演，该方法节约了内存空间，同时容易实现大规模并行运算，因此能够适用于实际的生产工作。王毓玮等提出了一种基于包络的多尺度反演方法，该方法利用包络反演为混合域弹性波全波形反演提供含有大尺度信息的初始模型，大大缓解了混合域弹性波全波形反演的“周波跳跃”现象。

基于包络的多尺度反演方法中首先是通过求取地震记录包络，并利用包络来进行反演得到大尺度的初始模型，然后利用得到的初始模型，应用Sirgue的混合域全波形反演方法实现速度模型的高精度建立。包络反演是建立在地震记录全频带基础上的，如果要进行分频段包络反演需要对地震记录做额外的低通滤波处理，这过程中的低通滤波器的设置同样会影响包络反演的效果。因此，这种方案虽然在一定程度上降低了弹性波全波形反演对低频信息的依赖，但是额外引入的低通滤波器却增加了工作量，使多尺度反演的分频框架不自然。

发明内容

本发明为解决上述现有技术方案中基于包络的弹性波全波形反演技术中分频反演需要额外引入低通滤波器的问题，提供了一种基于自然分频时频域能量谱的弹性波反演方法，减少了误差的产生，降低了反演陷入局部极值的概率，提高了速度模型的精度。

本发明采用的技术方案是：一种基于自然分频时频域能量谱的弹性波反演方法，包括以下步骤：

S1：输入实测三分量地震记录、从实测数据的道头中读入野外观测系统、读入地震子波wavelet、反演参数和初始纵波速度模型mvp1(0)、横波层速度模型mvs1(0)，依据反演参数对计算空间进行网格剖分；

S2：利用设置的参数和当前模型进行进行弹性波方程正演模拟，得到当前速度模型下的x和z分量的模拟地震记录currecordx、currecordz；

S3：读入x和z分量的实测地震记录orirecordx、orirecordz，确定时频域反演的频率范围[ω_min,ω_max]、设置反演频率间隔ω_gap，反演频率点个数nω，对x和z分量模拟记录做单频S变换得到curx、curz，对x和z分量实测地震记录做单频S变换得到orix、oriz，分别取curx、curz、orix、oriz的振幅的平方，得到各自的单频能量谱；

S4：利用自然分频的时频域能量谱弹性波全波形反演进行多尺度反演获得含有大尺度信息的纵横波层速度模型；

S5：将含有大尺度信息的速度模型作为初始模型进行混合域多尺度弹性波全波形反演得到最终的高精度反演结果。

本技术方案中，通过S变换将地震记录变换到时频域并利用不同单频点的能量谱进行反演，省去了常规包络反演中设计低通滤波器的步骤，减少了工作量，降低了工作难度，同时利用混合域多尺度全波形反演从低频段逐步反演到高频段，降低了反演陷入局部极值的概率，提高了速度模型的精度。

优选地，步骤S4利用自然分频的时频域能量谱弹性波全波形反演进行多尺度反演获得含有大尺度信息的纵横波层速度模型，具体步骤如下：

假设当前频率点号f，当前频率点下的最大迭代次数为iter1，第f个频率点的第i次迭代的初始纵波模型mvp1(i,f)、初始横波模型mvs1(i,f)，根据S3分别得到curx、curz、orix、oriz的单频能量谱|S_curx(t,f)|²、|S_curz(t,f)|²、|S_orix(t,f)|²、|S_oriz(t,f)|²，计算全波形反演的伴随震源adjx、adjz为

利用弹性波方程将adjx，adjz反传求得伴随波场，根据梯度类或牛顿类方法利用伴随波场与模拟波场求取梯度，确定此次迭代的纵波更新方向directionvp1,横波更新方向directionvs1,确定更新步长a；

纵波模型更新公式mvp1(i+1,f)＝mvp1(i,f)+a*directionvp1；

横波模型更新公式mvs1(i+1,f)＝mvs1(I,f)+a*directionvs1；

更新速度模型，判断频率点f的当前迭代次数i是否等于iter1，若等于，则进入下一个频率点，若不等于，则进行当前频率点的下一次迭代；判断频率点f是否等于nω，若等于，输出最终纵波速度模型和横波速度模型，则f＝f+1,i＝0,进入下一个频率点的迭代。

优选地，步骤S5包括将含有大尺度信息的速度模型作为初始模型进行混合域多尺度弹性波全波形反演获得高精度的速度模型，其具体步骤为：

S51：设置每一个频带范围[ifmin ifmax]，频带个数groupnum,频带中频率个数fnum,当前频带序号k,读入含有大尺度信息的初始模型minit,震源子波wavelet，设置横向网格间距dx,纵向网格间距dz,时间采样间隔dt，总时间采样点Tn,横向网格点数Xn,纵向网格点数Zn，最大迭代次数itermax，当前迭代数j，假设模型迭代过程中第k个频带的第j次的纵波模型、横波模型分别为mvp(k,j)、mvs(k,j)；

S52：利用wavelet、minit和设置的参数进行弹性波方程正演得到x、z分量的模拟记录CURX、CURZ并记录每一时刻的波场值的x分量和z分量ux、uz，读入实测记录orix、oriz；

S53：根据当前频带范围[ifmin ifmax]对ux、uz做离散傅里叶变换得到fux、fuz,模拟记录与实测记录分别做差得到记录残差Residualx、Residualz；

S54：根据当前频带范围[ifmin ifmax]对ux、uz做离散傅里叶变换得到fux、fuz,对CURX、orix做离散傅里叶变换并做差得到Residualx，对CURZ、oriz做离散傅里叶变换并做差得到Residualz；

S55：利用fux、fuz、frux、fruz互相关求取梯度gradvp、gradvs；

S56：梯度预处理确定更新方向directionvp、directionvs，计算更新步长a；

S57：更新模型：mvp(k,j+1)＝mvp(k,j)+a*directionvp；

mvs(k,j+1)＝mvs(k,j)+a*directionvs；

S58：判断j是否等于itermax；<1>如果是，进入步骤S59；<2>否则，令j＝j+1,进入步骤S52；

S59：判断k是否等于groupnum；<1>如果是，输出最终纵波速度模型和横波速度模型；<2>否则，令k＝k+1,j＝0,进入步骤S52。

与现有技术相比，有益效果是：本发明提出了基于时频域能量谱弹性波反演的初始模型构建方法，通过S变换将地震记录变换到时频域并取其最低可用频率的能量谱进行反演，省去了常规包络反演中设计低通滤波器的步骤，减少了误差的产生，同时利用混合域多尺度全波形反演从低频段逐步反演到高频段，降低了反演陷入局部极值的概率，提高了速度模型的精度。

附图说明

图1是本发明一种基于自然分频时频域能量谱的弹性波反演方法的流程图；

图2是本发明一种基于自然分频时频域能量谱的弹性波反演方法的详细框架流程图；

图3是本发明一种基于自然分频时频域能量谱的弹性波反演方法的精确速度模型图；

图4是本发明一种基于自然分频时频域能量谱的弹性波反演方法的初始线性模型图；

图5是本发明一种基于自然分频时频域能量谱的弹性波反演方法的传统混合域多尺度全波形反演结果图；

图6是本发明一种基于自然分频时频域能量谱的弹性波反演方法的时频域能量谱反演结果图；

图7是本发明一种基于自然分频时频域能量谱的弹性波反演方法的能量谱基础是进一步进行混合域多尺度反演效果图。

具体实施方式

附图仅用于示例性说明，不能理解为对本专利的限制；为了更好说明本实施例，附图某些部品会有省略、放大或缩小，并不代表实际产品的尺寸；对于本领域技术人员来说，附图中某些公知结构及其说明可能省略是可以理解的。附图中描述位置关系仅用于示例性说明，不能理解为对本专利的限制。

本发明实施例的附图中相同或相似的标号对应相同或相似的部品；在本发明的描述中，需要理解的是，若有术语“上”、“下”、“左”、“右”“长”“短”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此附图中描述位置关系的用语仅用于示例性说明，不能理解为对本专利的限制，对于本领域的普通技术人员而言，可以根据具体情况理解上述术语的具体含义。

下面通过具体实施例，并结合附图，对本发明的技术方案作进一步的具体描述：

实施例1

如图1-图2所示为一种基于自然分频时频域能量谱的弹性波反演方法的实施例1，包括以下步骤：

S1：输入实测三分量地震记录、从实测数据的道头中读入野外观测系统、读入地震子波wavelet、反演参数和初始纵波速度模型mvp1(0)、横波层速度模型mvs1(0)，依据反演参数对计算空间进行网格剖分；

S2：利用设置的参数和当前模型进行进行弹性波方程正演模拟，得到当前速度模型下的x和z分量的模拟地震记录currecordx、currecordz；

S3：读入x和z分量的实测地震记录orirecordx、orirecordz，确定时频域反演的频率范围[ω_min,ω_max]、设置反演频率间隔ω_gap，反演频率点个数nω，对x和z分量模拟记录做单频S变换得到curx、curz，对x和z分量实测地震记录做单频S变换得到orix、oriz，分别取curx、curz、orix、oriz的振幅的平方，得到各自的单频能量谱；

S4：利用自然分频的时频域能量谱弹性波全波形反演进行多尺度反演获得含有大尺度信息的纵横波层速度模型；

S5：将含有大尺度信息的速度模型作为初始模型进行混合域多尺度弹性波全波形反演得到最终的高精度反演结果。

本实施例的有益效果：本技术方案中，通过S变换将地震记录变换到时频域并利用不同单频点的能量谱进行反演，省去了常规包络反演中设计低通滤波器的步骤，减少了工作量，降低了工作难度，同时利用混合域多尺度全波形反演从低频段逐步反演到高频段，降低了反演陷入局部极值的概率，提高了速度模型的精度。

实施例2

一种基于自然分频时频域能量谱的弹性波反演方法的实施例2，在实施例1的基础上，对实施例1中的步骤进一步限定。

具体的，步骤S4利用自然分频的时频域能量谱弹性波全波形反演进行多尺度反演获得含有大尺度信息的纵横波层速度模型，具体步骤如下：

假设当前频率点号f，当前频率点下的最大迭代次数为iter1，第f个频率点的第i次迭代的初始纵波模型mvp1(i,f)、初始横波模型mvs1(i,f)，根据S3分别得到curx、curz、orix、oriz的单频能量谱|S_curx(t,f)|²、|S_curz(t,f)|²、|S_orix(t,f)|²、|S_oriz(t,f)|²，计算全波形反演的伴随震源adjx、adjz为

利用弹性波方程将adjx，adjz反传求得伴随波场，根据梯度类或牛顿类方法利用伴随波场与模拟波场求取梯度，确定此次迭代的纵波更新方向directionvp1,横波更新方向directionvs1,确定更新步长a；

纵波模型更新公式mvp1(i+1,f)＝mvp1(i,f)+a*directionvp1；

横波模型更新公式mvs1(i+1,f)＝mvs1(I,f)+a*directionvs1；

更新速度模型，判断频率点f的当前迭代次数i是否等于iter1，若等于，则进入下一个频率点，若不等于，则进行当前频率点的下一次迭代；判断频率点f是否等于nω，若等于，输出最终纵波速度模型和横波速度模型，则f＝f+1,i＝0,进入下一个频率点的迭代。

具体的，步骤S5包括将含有大尺度信息的速度模型作为初始模型进行混合域多尺度弹性波全波形反演获得高精度的速度模型，其具体步骤为：

S51：设置每一个频带范围[ifmin ifmax]，频带个数groupnum,频带中频率个数fnum,当前频带序号k,读入含有大尺度信息的初始模型minit,震源子波wavelet，设置横向网格间距dx,纵向网格间距dz,时间采样间隔dt，总时间采样点Tn,横向网格点数Xn,纵向网格点数Zn，最大迭代次数itermax，当前迭代数j，假设模型迭代过程中第k个频带的第j次的纵波模型、横波模型分别为mvp(k,j)、mvs(k,j)；

S52：利用wavelet、minit和设置的参数进行弹性波方程正演得到x、z分量的模拟记录CURX、CURZ并记录每一时刻的波场值的x分量和z分量ux、uz，读入实测记录orix、oriz；

S53：根据当前频带范围[ifmin ifmax]对ux、uz做离散傅里叶变换得到fux、fuz,模拟记录与实测记录分别做差得到记录残差Residualx、Residualz；

S54：根据当前频带范围[ifmin ifmax]对ux、uz做离散傅里叶变换得到fux、fuz,对CURX、orix做离散傅里叶变换并做差得到Residualx，对CURZ、oriz做离散傅里叶变换并做差得到Residualz；

S55：利用fux、fuz、frux、fruz互相关求取梯度gradvp、gradvs；

S56：梯度预处理确定更新方向directionvp、directionvs，计算更新步长a；

S57：更新模型：mvp(k,j+1)＝mvp(k,j)+a*directionvp；

mvs(k,j+1)＝mvs(k,j)+a*directionvs；

S58：判断j是否等于itermax；

如果是，进入步骤S59；如果否，令j＝j+1,进入步骤S52；

S59：判断k是否等于groupnum；

如果是，输出最终纵波速度模型和横波速度模型；如果否，令k＝k+1,j＝0,进入步骤S52。

本实施例的有益效果：本发明提出了基于时频域能量谱弹性波反演的初始模型构建方法，通过S变换将地震记录变换到时频域并取其最低可用频率的能量谱进行反演，省去了常规包络反演中设计低通滤波器的步骤，减少了误差的产生，同时利用混合域多尺度全波形反演从低频段逐步反演到高频段，降低了反演陷入局部极值的概率，提高了速度模型的精度。

实施例3

一种基于自然分频时频域能量谱的弹性波反演方法的实施例3，如图3-图7所示，在实施例1或实施例2的基础上，举例说明本方法的优势。

具体的，如图3所示，展示了所用准确模型，采用8Hz主频的Ricker子波激发，横向网格步长和纵向网格步长均为10m，时间采样间隔为1ms，总采样时间为3s,炮间隔100m,全排列接收，炮点和检波点位于地表，得到地震记录，滤去地震记录中的2.5Hz以下的低频信息得到低频缺失的实测记录，注意纵波速度与横波速度之比固定为1.6。

如图4所示，展示了所用的1D线性初始模型，可以看出在模型中不包含任何的构造信息，这对弹性波全波形反演来说是一个很大的挑战。

如图5所示，展示了直接利用图4进行混合域多尺度弹性波全波形反演的结果，反演过程中利用L_BFGS算法进行优化加速，可以看出得益于多尺度方法的运用，模型浅部构造得到了一定的恢复，但是由于缺失低频并且初始模型较差，常规混合域多尺度弹性波全波形反演所受周期跳跃现象的影响严重，因此在模型深部以及箭头所指的断裂面附近反演效果较差，而由于横波速度模型反演对于低频信息和初始模型的依赖性更强，因此在缺失低频且初始模型较差时，其反演效果更差。

如图6所示，展示了时频域能量谱反演得到的含有大尺度构造的速度模型，在进行时频域能量谱反演实验时，利用L_BFGS算法进行优化加速，可以看到相比于图4，采用本发明基于时频域能量谱的多尺度弹性波全波形反演提供的初始模型更加精确，在低频缺失情况下恢复了模型的大尺度构造，这为后续的混合域多尺度反演提供了良好的初始模型，并且本发明在反演过程中并没有额外引入低通滤波器，减少了工作量，降低了工作难度。

如图7所示，展示了基于自然分频时频域能量谱反演的结果进一步开展混合域多尺度弹性波全波形反演的结果，可以看到，模型中的断裂面、高速层等都得到有效的重建。相比图5，图7在相同区域内的反演精度更高，分辨率也更高，并且可以看出本发明对模型深部的反演效果提升尤为明显。这主要得益于本发明所采用的自然分频时频域多尺度能量谱反演提供的较好的初始模型，这在一定程度上避免了周期跳跃现象。本发明对低频缺失下的纵波速度模型和横波速度模型的反演效果均有提升，并且对横波反演效果的提升尤为明显。

显然，本发明的上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明权利要求的保护范围之内。

Claims (2)

1.一种基于自然分频时频域能量谱的弹性波反演方法，其特征在于，包括以下步骤：

S1：输入实测三分量地震记录、从实测数据的道头中读入野外观测系统、读入地震子波wavelet、反演参数和初始纵波速度模型mvp1(0)、横波层速度模型mvs1(0)，依据反演参数对计算空间进行网格剖分；

S2：利用设置的参数和当前模型进行进行弹性波方程正演模拟，得到当前速度模型下的x和z分量的模拟地震记录currecordx、currecordz；

S3：读入x和z分量的实测地震记录orirecordx、orirecordz，确定时频域反演的频率范围[ω_min,ω_max]、设置反演频率间隔ω_gap，反演频率点个数nω，对x和z分量模拟记录做单频S变换得到curx、curz，对x和z分量实测地震记录做单频S变换得到orix、oriz，分别取curx、curz、orix、oriz的振幅的平方，得到各自的单频能量谱；

S4：利用自然分频的时频域能量谱弹性波全波形反演进行多尺度反演获得含有大尺度信息的纵横波层速度模型；

S5：将含有大尺度信息的速度模型作为初始模型进行混合域多尺度弹性波全波形反演得到最终的高精度反演结果。

2.根据权利要求1所述的一种基于自然分频时频域能量谱的弹性波反演方法，其特征在于，在所述S4中，具体实施方式为：

假设当前频率点号f，当前频率点下的最大迭代次数为iter1，第f个频率点的第i次迭代的初始纵波模型mvp1(i,f)、初始横波模型mvs1(i,f)，根据S3分别得到curx、curz、orix、oriz的单频能量谱|S_curx(t,f)|²、|S_curz(t,f)|²、|S_orix(t,f)|²、|S_oriz(t,f)|²，计算全波形反演的伴随震源adjx、adjz为

利用弹性波方程将adjx，adjz反传求得伴随波场，根据梯度类或牛顿类方法利用伴随波场与模拟波场求取梯度，确定此次迭代的纵波更新方向directionvp1,横波更新方向directionvs1,确定更新步长a；

纵波模型更新公式mvp1(i+1,f)＝mvp1(i,f)+a*directionvp1；

横波模型更新公式mvs1(i+1,f)＝mvs1(I,f)+a*directionvs1；

更新速度模型，判断频率点f的当前迭代次数i是否等于iter1，若等于，则进入下一个频率点，若不等于，则进行当前频率点的下一次迭代；判断频率点f是否等于nω，若等于，输出最终纵波速度模型和横波速度模型，则f＝f+1,i＝0,进入下一个频率点的迭代。