一种基于声波方程的非线性优化隐式时空域有限差分数值模拟方法

摘要

本发明公开了一种基于声波方程的非线性优化隐式时空域有限差分数值模拟方法，包括以下步骤：(1)读取参数；(2)求解二阶时间导数的高阶离散差分格式，并计算高阶差分系数；(3)采用空间隐式离散格式求解二阶空间导数，求解隐式差分系数；(4)构建同时具有时间高阶和空间隐式的差分递推格式，得到时空域频散关系；(5)采用非线性优化算法对时空域差分系数进行求解；(6)采用混合吸收边界条件对反射进行吸收，按照波动方程进行递推，得到任意时刻的波场及整个地震记录；(7)记录波场快照，输出地震记录并结束。本发明提出的方法在保证精度的前提下可以显著减小计算时间，为逆时偏移和全波形反演方法高效地提供精度更高的波场。

说明书

技术领域

本发明属于地球物理勘探数据处理技术领域，涉及一种基于求解声波方程的高精度、高效率非线性优化隐式时空域有限差分数值模拟方法。

背景技术

目前基于声波波动方程的逆时偏移方法相比基于单程波方程的偏移成像方法和基于射线理论的克西霍夫成像方法，在处理盐丘等复杂构造方面已经显示了巨大的优势。基于声波方程的全波形反演方法在反演储层参数方面也越来越受到地球物理学家的重视。其主要原因是这两种方法充分利用了地震波的运动学和动力学特征，在成像和反演中考虑了一次波、多次波、反射波和折射波等各种信息，能够最大程度的揭示地下构造信息和储层参数。全波形反演和逆时偏移方法都涉及波动方程的求解，计算时间长，计算量巨大，某种程度上制约着这两种技术的发展。因此高精度、高效率的数值模拟方法就显得尤其重要。有限差分法相比有限元法、伪谱法等其他数值模拟方法，具有操作简单、计算效率高、易于并行和所需内存小等诸多优势，被广泛应用于声波方程的数值求解。

目前的有限差分法主要缺点是易于频散且稳定性差。频散按照离散对象可以分为时间差分频散和空间差分频散。时间频散可以采用较小的时间步长来压制，但这同样带来巨大的计算时间。而直接采用高阶差分近似时间导数会导致传播过程不稳定。空间频散可以通过增大算子长度或减小空间步长来压制。但这两种策略都会增大计算时间。空间频散还可以通过采用隐式差分法进行压制。相比显式差分法，隐式差分法可以采用较小的算子长度来达到与显式差分一致的精度，因此有助于减小计算时间。然而，由于时间频散的存在，空间精度的提高并不会有效压制整体频散。

发明内容

本发明的目的是为减小传统隐式有限差分方法时间频散严重、稳定性差的问题而提供的一种基于声波方程的非线性优化隐式时空域有限差分数值模拟方法，本发明提出的方法在保证精度的前提下可以显著减小计算时间，为逆时偏移和全波形反演方法高效地提供精度更高的波场。

本发明的一种基于声波方程的非线性优化隐式时空域有限差分数值模拟方法，包括以下几个步骤：

(1)读取参数，所述参数包括正演模拟所需要的速度模型参数文件、有限差分算子长度、非线性迭代求解差分系数时需要的最大允许误差、子波函数及子波主频、正演所采用的时间与空间步长及地震记录时长；

(2)求解二阶时间导数的高阶离散差分格式，并计算高阶差分系数；

(3)借鉴Claerbout(1985)提出的求解二阶空间导数隐式差分格式，在常规显式高阶差分的基础上，分母中引入二阶中心差分格式，构造空间隐式差分格式求解二阶空间导数，用于减小有限差分算子长度；为进一步减小频散，基于空间频散关系，求解隐式差分系数；

(4)基于所述时间高阶离散格式和空间隐式离散格式，将差分格式带入声波方程，得到同时具有时间高阶和空间隐式的差分递推格式；利用平面波分析理论，将平面波方程带入频散关系，化简得到时空域频散关系；

(5)以时间和空间差分离散格式中的差分系数为变量，通过极小化方程频散关系，建立关于时间和空间差分系数的非线性目标函数；以步骤(2)和步骤(3)计算得到的差分系数为初始值，采用Solvopt非线性优化算法对时间及空间差分系数进行迭代求解；

(6)采用混合吸收边界条件对边界反射进行吸收，按照波动方程的离散格式进行递推，得到任意时刻的波场及整个地震记录；

(7)记录波场快照，输出地震记录并结束。

步骤(2)中，求解二阶时间导数时采用中心差分，形式如下：

其中，为压力场，τ为时间步长，t为时间；为了进一步提高时间离散精度，将菱形差分算子引入到公式(1)中，得出如下的时间高阶离散公式：

其中，h为空间采样步长，v为速度，a_m,n为有限差分系数；N为菱形算子中的差分算子长度，其值越大，时间差分近似精度越高，但增加的计算量也增大；n和m为菱形算子中的坐标位置序列。

上述有限差分系数a_m,n通过对公式(2)采用泰勒展开法计算得到，公式如下：

其中，j＝2,3,…,N，ξ＝1,2,…,int(j/2)，r为库朗数，r＝vτ/h，int是取整运算；通过采用公式(2)对时间导数进行离散，并采用公式(3)所示的差分系数，时间离散精度达到高阶，用于压制时间频散。

步骤(3)中，采用如下空间隐式离散格式对二阶空间导数进行近似，形式如下：

其中，a′₀、a′_m和b为隐式差分系数，h为空间步长，x为方向坐标，M为差分算子长度；公式(4)对应的空间频散关系如下：

其中，β＝k_xh，k_x为水平方向波数；极小化公式(4)对应的空间频散关系，目标函数为关于差分系数的线性凸函数，最优差分系数通过求解如下公式得到：

其中，n＝1,2,…,M+1，d＝[d₁,d₂,…,d_M,d_M+1]^T＝[a′₁,a′₂,…,a′_M-1,a′_M,b]^T，q为介于0和1之间的自然数，控制着积分范围，

步骤(4)中，在求解时间导数时，采用公式(2)计算，求解空间导数时，采用公式(4)计算，将公式(2)和公式(4)带入到声波方程并进行平面波分析，过程如下：

均匀模型中，声波方程压力P满足如下平面波方程，

其中，ω是角频率，是虚数单位，k_z为垂直方向波数；m、n和l为位置和时间序列；

将公式(8)带入公式(2)并进行化简得到如下公式：

同样，对空间导数公式进行分析，得到如下频散关系：

将上述(9)、(10)和(11)带入声波方程并简化和重组得到如下时空域频散关系，

其中，ω为角频率，k_x和k_z分别为沿x和z方向的波数。

步骤(5)具体的步骤如下：

首先，差分系数a′_m和b由公式(6)计算得到，差分系数a_m,0和a_m,n由公式(3)计算得到；基于时空域频散关系(12)，以上述求得的差分系数为初始输入，利用非线性优化算法求解时空域差分系数：

定义如下目标函数，

其中，θ为传播角度，q为介于0和1间的实数，

公式(13)中的解通过非线性优化算法计算得到，公式如下：

(a_m,0,a_m,n,a′_m,b,r)_optimal＝Minimize[J₂(a_m,0,a_m,n,a′_m,b,r)](16)

上述非线性优化算法具体采用的是SolvOpt非线性优化算法求取，采用如下策略来减小计算差分系数所需要的运算时间：

假定模拟对象的速度模型(速度模型是正演模拟的对象，作为输入输入到程序中)参数中的速度取值范围为[v₁,v₂],则在计算差分系数时，每间隔dv计算差分系数，形式如下：

位于[v₁,v₁+dv]之间的速度对应的差分系数与v₁处的差分系数一致。尽管这样操作一定程度上损失了计算精度，但由于该方法大大降低了计算时间，可以显著提高计算效率。且复杂模型的计算发现，这种策略引起的计算误差并不大。

步骤(6)中，整个计算区域被划分为内部区域、边界区域和过渡区域，具体的方法如下：

(6-1)利用离散公式(2)和(4)离散时间及空间导数，采用公式(17)计算差分系数，进行波场递推，得到所有区域处的任意时刻的双程波波场；

(6-2)求解声波单程波方程，计算边界区域和过渡区域处的单程波方程，计算公式及方法如下：

以左边界和左上角为例，左边界单程波方程为：

其中，x和z代表坐标；左上角单程波方程为：

其他边界处及角点处单程波类似得到，并进行对应求解；采用上述吸收边界条件并按照上述方法进行递推，得到任意时刻波场值；

(6-3)计算最终波场；

边界处波场为单程波场值，内部区域处波场为双程波场值，中间过渡区域处波场通过如下线性组合得到：

P＝B_iP^two+(1-B_i)P^one,(20)

其中，P^two为双程波场，P^one为单程波场，B_i为权重系数，B_i＝(i-1)/N_b，i＝2,…,N_b，N_b为过渡区域吸收边界层数；过渡区域从内到外，B_i由(N_b-1)/N_b变化到1/N_b，如果N_b＝1，则此时混合吸收边界条件退化为单程波吸收边界条件。

本发明与现有技术相比其显著优点在于：

a)时间精度高，频散小。

b)稳定性好，模拟中可以采用更大的时间步长。

c)方法空间频散相比传统泰勒展开法被更好的压制，可以采用更小的算子长度，有效节省计算时间，计算效率高。

附图说明

图1为本发明的一种基于声波方程的非线性优化隐式时空域有限差分数值模拟方法工作流程示意图；

图2(a)为常规隐式差分方法频散曲线随时间步长变化图；

图2(b)为本发明方法频散曲线随时间步长变化图；

图3为本发明方法稳定性条件曲线与常规方法对比图；

图4(a)为常规隐式差分法计算得到的均匀模型(1200m,2250m)处震动曲线同参考解对比图；

图4(b)为本方法计算得到的均匀模型(1200m,2250m)处震动曲线同参考解对比图；

图4(c)为常规隐式差分法及本发明方法模拟震动曲线同参考解之间的统计均方根误差对比图；

图5为SEG/EAGE盐丘速度模型；

图6(a)为SEG/EAGE盐丘模型2.4s时刻参考波场快照；

图6(b)为SEG/EAGE盐丘模型应用常规隐式差分方法后2.4s时刻波场快照；

图6(c)为SEG/EAGE盐丘模型应用本发明方法后2.4s时刻波场快照；

图7(a)为SEG/EAGE盐丘模型应用常规隐式差分方法后(500m,4500m)处地震记录震动图与参考解对比图；

图7(b)为SEG/EAGE盐丘模型应用本发明方法后(500m,4500m)处地震记录震动图与参考解对比图。

具体实施方式

下面结合附图和实施例对本发明的具体实施方式作进一步的详细说明。

本发明提供了一种基于声波方程的非线性优化隐式时空域有限差分数值模拟方法，该方法首先在常规时间二阶离散格式的基础上，加入了菱形差分算子，推导了具有高阶精度的时间离散格式及差分系数。求解空间导数时，本发明采用优化的隐式差分方法，在保证精度的前提下减小算子长度。结合高阶隐式空间差分方法，本发明得出了同时具有高阶时间精度和隐式空间精度的递推格式及其时空域频散关系。从时空域频散关系出发，本发明采用非线性优化算法，提出了计算时空域差分系数的方法和策略。相比常规方法，本发明方法精度高，稳定性好。为避免在复杂模型中计算每个速度处的差分系数，本发明进一步提供了一种节省计算时间的选择速度策略。

如图1所示，本发明提出的一种基于声波方程的非线性优化隐式时空域有限差分数值模拟方法包括如下具体步骤：

(1)读取参数，包括正演模拟所需要的速度模型参数文件、有限差分算子长度(时间导数和空间导数)、非线性迭代求解差分系数时需要的终止条件、子波函数及子波主频、正演所采用的时间和空间步长及地震记录时长等参数；

(2)针对时间导数的高阶离散差分格式构造和差分系数求解：

传统方法在求解2阶时间导数时采用中心差分，形式如下：

其中，为压力场，τ为时间步长，t为时间。为了进一步提高时间离散精度，本发明将菱形差分算子引入到公式(1)中，提出了如下的时间高阶离散公式：

其中，h为空间采样步长，v为速度，N为菱形算子中的差分算子长度，其值越大，时间差分近似精度越高，但增加的计算量也增大。n和m为菱形算子中的坐标位置序列。有限差分系数a_m,n可以通过对公式(2)采用泰勒展开法得到，公式如下：

其中，r为库朗数，r＝vτ/h，int是取整运算。通过采用公式(2)对时间导数进行离散，并采用公式(3)所示的差分系数，时间离散精度可以达到高阶，有效压制时间频散。

(3)利用隐式差分方法求解空间导数，为进一步减小频散，基于空间频散关系，通过优化方法求解隐式差分系数：

针对二阶空间导数，为有效减小算子长度，本发明采用如下隐式差分离散格式对空间二阶导数进行近似，形式如下：

其中，a′₀、a′_m和b为隐式差分系数。h为空间步长，x为方向坐标，M为差分算子长度。常规方法在求取差分系数a′₀、a′_m和b时采用泰勒展开法，精度低，易于频散。本文采用优化方法求解，公式(4)对应的空间频散关系如下：

其中，β＝k_xh，k_x为水平方向波数。极小化公式(4)对应的空间频散关系，目标函数为关于差分系数的线性凸函数，可以通过求解如下公式得到：

其中，d＝[d₁,d₂,…,d_M,d_M+1]^T＝[a′₁,a′₂,…,a′_M-1,a′_M,b]^T，q为介于0和1之间的自然数，控制着积分范围，

(4)构建同时具有时间高阶和空间隐式的差分递推格式，并得到时空域频散关系：

在求解时间导数时，采用公式(2)计算，求解空间导数时，采用公式(4)计算，将公式(2)和公式(4)带入到声波方程并进行平面波分析方法，可以得到如下时空域频散关系，

其中，ω为角频率，k_x和k_z分别为沿x和z方向的波数。

(5)从时空域频散关系出发，以上述方法求得的差分系数为初始猜测，采用非线性优化算法对时空域差分系数进行迭代优化求解：

首先，差分系数a′_m和b由公式(6)计算得到，差分系数a_m,0和a_m,n由公式(3)计算得到。在通过上述方法得到差分系数a_m,0、a_m,n、a′_m和b后，频散误差仍然可以通过非线性优化算法进行进一步降低，定义如下目标函数，

其中，θ为传播角度，q为介于0和1间的实数，

公式(9)可以通过非线性优化算法计算得到，公式如下：

(a_m,0,a_m,n,a′_m,b,r)_optimal＝Minimize[J₂(a_m,0,a_m,n,a′_m,b,r)].(12)

非线性优化算法具体采用的是SolvOpt非线性优化算法求取，SolvOpt是一种非线性局部优化算法，目标函数收敛速度快，有助于快速求解非线性优化差分系数。

对于复杂模型，上述运算需要对每个速度值进行计算，由于涉及迭代求解，运算时间很长。本发明提出采用如下策略来减小计算差分系数所需要的运算时间。

假定模拟对象的速度模型中涉及的速度取值范围为[v₁,v₂],则在计算差分系数时，每间隔dv计算差分系数，形式如下：

位于[v₁,v₁+dv]之间的速度对应的差分系数与v₁处的差分系数一致，其他速度区间与该区间计算方法类似。通过上述策略，计算差分系数所需要的时间大大降低，同时用于存储差分系数的内存也相应降低。同时，复杂模型表明，该策略并不会显著降低差分方法的精度。

(6)采用适当的吸收边界条件，进行波场延拓：

求得差分系数后，按照波动方程递推格式进行递推，得到任意时刻的波场及整个地震记录。本发明采用Liu和Sen(2010)提出的混合吸收边界条件对反射进行吸收。

(7)记录波场快照，输出地震记录并终止程序。

以下为本发明的一个实施例，说明基于一种基于声波方程的非线性优化隐式时空域有限差分数值模拟方法的实现过程和最终结果。为了便于阐述，在下文中，将传统方法记做ISD，将本发明提出的方法记做ITSD-2。

图2(a)和图2(b)为本发明提出的方法与传统方法频散误差对比分析，图中纵坐标为数值模拟速度与真实速度比值，反映了方法的精度。当纵坐标为1.0时，此时方法无频散。纵坐标偏离1.0的位置越远，则频散越大。图中，速度参数为1500m/s，空间采样步长为15m，算子长度M＝6。图中显示的是传播方向为45度时采用不同时间步长时频散曲线。对比发现，相比传统方法，本发明中提出的方法显著改进了传统方法的精度，压制了数值频散。

图3为本发明提出的方法与传统方法稳定性对比分析图，图中纵坐标为稳定性因子，当库朗数大于该稳定性因子时，差分方法不稳定。该值越大，方法越稳定。对比发现，传统方法稳定性最差，稳定性因子最小。本发明提出的方法稳定性更好，有助于方法采用更大的时间步长。

图4a和图4b为均匀模型(1200m,2250m)处模拟波场震动图与参考解对比图。模型大小为4500m×4500m，空间步长为15m，时间步长为3ms，采用的子波为主频是18Hz的雷克子波，在模型中心激发。参考解为由Cagniard-de Hoop(De Hoop,1960)方法计算得到。本发明中方法采用M＝6和N＝2。对比发现，传统方法频散误差明显，与参考解吻合差。本发明提出的方法则显著提高了精度，与参考解吻合更好。图4c为各方法与参考解之间的统计均方根误差，明显本发明提出的方法误差更小。

图5为地震勘探数值模拟中常用到的SEG/EAGE盐丘模型。模型大小为5000m×15000m，空间步长为25m，采用的子波为主频是15Hz的雷克子波，在地表中心激发，时间步长为2.5ms。本发明中方法采用M＝8和N＝2。图6(a)至图6(c)为各方法在2.4s时刻时传播波场快照，参考解为常规显式时空域方法采用很小的时间步长和很大的算子长度计算得到。同参考解对比发现，本发明中方法明显具有更高的精度，白色箭头所指区域频散相比传统方法明显小。图7(a)和图7(b)为各方法在(500m,4500m)处地震记录震动图相比参考解对比图，图中显示的误差值为统计的均方根误差，明显本发明提出的方法均方根误差更小，且非线性优化方法相比线性优化方法精度更高，误差更小。

为降低频散误差，传统隐式差分方法采用1ms时间步长时计算的均方根误差为1.1824e^-4，消耗的CPU时间为403.58s。本发明提出的方法当采用2.5ms时间步长时，计算的均方根误差为1.0206e^-4，消耗的CPU时间为240.27s。可以发现，本发明提出的方法在保持了精度的前提下，显著降低了计算时间，效率得到大大提高。

本方法可以在保证精度的前提下，显著提高计算效率，为计算量巨大的逆时偏移方法和全波形反演方法提供便利。本发明的具体实施方式中未涉及的说明属于本领域公知的技术，可参考公知技术加以实施。

以上具体实施方式及实施例是对本发明提出的一种基于声波方程的非线性优化隐式时空域有限差分数值模拟方法技术思想的具体支持，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在本技术方案基础上所做的任何等同变化或等效的改动，均仍属于本发明技术方案保护的范围。