本发明属于地震波场数值模拟技术领域，涉及一种基于求解声波方程的高精度、高效率线性优化隐式时空域有限差分数值模拟方法。
背景技术：
目前逆时偏移方法和全波形反演技术在处理如盐丘等复杂构造的成像和储层参数反演方面已经显示出了巨大的潜力，越来越受到地球物理学家的重视。但逆时偏移和全波形反演方法都涉及波动方程的数值求解，在模型偏大尤其三维时，计算量大，耗时长，某种程度上制约着这两种技术的发展。因此高精度、高效率的数值模拟方法就显得尤其重要。
有限差分法相比其他数值模拟方法，具有操作简单、计算效率高、易于并行和所需内存小等诸多优势，在逆时偏移和全波形反演中被广泛应用。但目前有限差分法存在频散严重、精度低和稳定性差等缺点。频散按照离散对象可以分为时间差分频散和空间差分频散。
技术实现要素：
本发明的目的是为克服传统隐式空间有限差分方法时间精度低、稳定性差和易于频散等问题而提供一种基于声波方程的线性优化隐式时空域有限差分数值模拟方法，具有更高精度和更好稳定性，本发明在保持空间隐式差分精度的前提下，提出了显著压制时间频散的线性优化时空域策略(即权利请求书中的第5步骤)，可以为逆时偏移和全波形反演方法提供频散更小、精度更高的波场，同时可以帮助提高声波方程数值模拟的精度和效率。
本发明的一种基于声波方程的线性优化隐式时空域有限差分数值模拟方法，包括以下几个步骤：
(1)读取参数，所述参数包括正演模拟所需要的速度模型参数文件、有限差分算子长度、子波函数及子波主频、正演所采用的时间与空间步长及地震记录时长；
(2)基于菱形差分算子和二阶时间中心差分，得到时间导数的时间高阶离散格式；
(3)借鉴claerbout(1985)提出的求解二阶空间导数隐式差分格式，在常规显式高阶差分的基础上，分母中引入二阶中心差分格式，构造空间隐式差分格式求解二阶空间导数，用于减小有限差分算子长度；进一步，基于最小二乘法，以l2范数目标函数极小化空间域频散关系，求解优化的隐式差分系数；
(4)基于所述时间高阶离散格式和空间隐式离散格式，带入声波方程，构建同时具有时间高阶和空间隐式的差分递推格式，得到时空域频散关系；
(5)固定空间差分系数，采用优化策略极小化整个递推格式频散关系，建立关于时间差分格式中差分系数的l2范数目标函数，利用线性优化算法对时间差分系数进行求解(详细内容见下文公式(11)～(16))；
(6)采用liu等(2010)提出的混合吸收边界条件对边界反射能量进行吸收，按照上述介绍的时间高阶和空间隐式差分格式离散波动方程并进行递推，得到任意时刻的波场及整个地震记录；
(7)记录波场快照，输出地震记录并结束。
步骤(2)中，求解二阶时间导数时采用中心差分，形式如下：
其中，
其中，h为空间采样步长，v为速度，am,n为有限差分系数；n为菱形算子中的差分算子长度，其值越大，时间差分近似精度越高，但增加的计算量也增大；n和m为菱形算子中的坐标位置序列。
步骤(3)中，采用如下空间隐式离散格式对二阶空间导数进行近似，形式如下：
其中，a′0、a′m和b为隐式差分系数，当b＝0时，公式(3)退化为显式差分格式；h为空间步长，x为方向坐标，m为差分算子长度；公式(3)对应的空间频散关系如下：
其中，β＝kxh，kx为水平方向波数；极小化公式(4)对应的空间频散关系，目标函数为关于差分系数的线性凸函数，最优差分系数通过求解如下公式得到：
其中，n＝1,2,…,m+1，d＝[d1,d2,…,dm,dm+1]t＝[a′1,a′2,…,a′m-1,a′m,b]t，q为介于0和1之间的自然数，控制着积分范围，
步骤(4)中，在求解时间导数时，采用公式(2)计算；求解空间导数时，采用公式(3)计算，将公式(2)和公式(3)带入到声波方程并进行平面波理论分析，过程如下：
均匀模型中，声波方程压力p满足如下平面波方程，
其中，ω是角频率，
将公式(7)带入公式(2)并进行化简得到如下公式：
同样，对空间导数公式进行分析，得到如下频散关系：
将上述(8)、(9)和(10)带入声波方程并简化和重组得到如下时空域频散关系，
其中，ω为角频率，kx和kz分别为沿x和z方向的波数，r为库朗数，r＝vτ/h。
步骤(5)中，首先，差分系数a′m和b由公式(5)计算得到，在计算差分系数am,0和am,n时保持固定，此时优化差分系数am,0和am,n是线性优化问题，直接线性求解，定义如下目标函数：
其中，d＝é[d1,…,dn+(n-1)n/2]t＝[a1,0,…,an,0,a1,1,…,a1,n-1,…,an-1,1]t，表示差分系数；θ为传播角度，q为介于0和1间的实数，波数kx和kz与波数k之间满足如下关系：
kx＝kcos(θ),kz＝ksin(θ),(13)ψ(β,θ,r)及
目标函数是关于差分系数am,0和am,n的线性函数，差分系数am,0和am,n通过求解如下公式得到，
步骤(6)中，整个计算区域被划分为内部区域、边界区域和过渡区域，具体的方法如下：
(6-1)利用公式(2)和公式(3)离散时间和空间导数，采用公式(5)和公式(16)计算差分系数，进行双程波方程求解，求解声波方程，计算所有区域处的双程波波场；
(6-2)求解单程波方程，计算边界区域和过渡区域处的单程波波场，计算公式及方法如下：
针对二阶声波方程，以左边界和左上角为例，左边界单程波方程为：
其中，x和z代表坐标；
左上角单程波方程为：
其他边界处及角点处单程波可以类似得到，并进行对应求解；采用上述吸收边界条件策略计算各位置处任意时刻波场值；
(6-3)计算最终波场：
边界处波场为单程波场值，内部区域处波场为双程波场值，中间过渡区域处波场通过如下线性组合得到：
p＝biptwo+(1-bi)pone,(19)
其中，ptwo为双程波场，pone为单程波场，bi为权重系数，bi＝(i-1)/nb，i＝2,…,nb，nb为过渡区域吸收边界层数，过渡区域从内到外，bi由(nb-1)/nb变化到1/nb；如果nb＝1，则此时混合吸收边界条件退化为单程波吸收边界条件。
本发明与现有技术相比其显著优点在于：
(a)时间精度高，频散小。在使用相同时间步长且步长偏大时，本发明方法可以得到精度更高的波场。
(b)稳定性好，模拟中可以采用更大的时间步长。
(c)时间精度提高的同时，空间隐式精度很好地保持。
附图说明
图1为本发明的一种基于声波方程的线性优化隐式时空域有限差分数值模拟方法工作流程示意图；
图2(a)为常规隐式空域方法频散曲线随时间步长变化图；
图2(b)为本发明方法当n＝2时频散曲线随时间步长变化图；
图2(c)为本发明方法当n＝3时频散曲线随时间步长变化图；
图3为本发明方法稳定性条件曲线与常规方法对比图；
图4(a)为均匀模型常规隐式差分方法模拟震动曲线同参考解对比图；
图4(b)为均匀模型本发明方法在n＝2时模拟震动曲线同参考解对比图；
图4(c)为本发明方法在n＝3时模拟震动曲线同参考解对比图；
图4(d)为本发明方法与常规隐式差分法模拟结果相对于参考解均方根误差对比图；
图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阶时间导数时采用中心差分，形式如下：
其中，
其中，h为空间采样步长，v为速度，am,n为有限差分系数。n为菱形算子中的差分算子长度，其值越大，时间差分近似精度越高，但增加的计算量也增大。n和m为菱形算子中的坐标位置序列。(3)针对空间导数，采用隐式差分格式求解以减小有限差分算子长度。为进一步减小频散，基于空间频散关系，通过优化方法求解隐式差分系数；
针对二阶空间导数，为有效减小算子长度，本发明采用如下空间隐式离散格式对空间二阶导数进行近似，形式如下：
其中，a′0、a′m和b为隐式差分系数，h为空间步长，x为方向坐标，m为差分算子长度。常规方法在求取差分系数a′m和b时采用泰勒展开法，精度低，易于频散。本文采用优化方法求解，公式(3)对应的空间频散关系如下：
其中，β＝kxh，kx为水平方向波数。极小化公式(4)左右两边的绝对误差，目标函数为关于差分系数的线性凸函数，可以通过求解如下公式得到：
其中，d＝[d1,d2,…,dm,dm+1]t＝[a′1,a′2,…,a′m-1,a′m,b]t，q为介于0和1之间的自然数，控制着积分范围，
(4)结合时间高阶离散格式和空间隐式离散格式，构建同时具有时间高阶和空间隐式的差分递推格式，得到时空域频散关系；
在求解时间导数时，采用公式(2)计算，求解空间导数时，采用公式(3)计算，将公式(2)和公式(3)带入到声波方程并进行平面波分析，可以得到如下时空域频散关系，
其中，ω为角频率，kx和kz分别为沿x和z方向的波数，r为库朗数，r＝vτ/h。
(5)从时空域频散关系出发，采用线性优化算法对时空域差分系数进行求解；
不同于显式差分方法，隐式差分方法频散关系形式更为复杂，直接同时优化am,0、am,n、a′m和b是强非线性问题。
首先，差分系数a′m和b由公式(5)计算得到，在计算差分系数am,0和am,n时保持固定，此时优化差分系数am,0和am,n是线性优化问题，可以直接线性求解，定义如下目标函数：
其中，d＝[d1,…,dn+(n-1)n/2]t＝[a1,0,…,an,0,a1,1,…,a1,n-1,…,an-1,1]t，表示差分系数，θ为传播角度，q为介于0和1间的实数，
目标函数(8)是关于差分系数am,0和am,n的线性函数，差分系数am,0和am,n可以通过求解如下公式得到，
(6)利用求取的差分系数，采用适当的吸收边界条件，进行波场延拓；
求得差分系数后，按照波动方程进行递推，得到任意时刻的波场及整个地震记录。波场延拓中必须要考虑的一个问题是边界反射。本发明采用liu和sen(2010)提出的混合吸收边界条件对边界反射进行吸收。
(7)记录波场快照，输出地震记录并终止程序。
以下为本发明的一个实施例，说明基于一种基于声波方程的线性优化隐式时空域有限差分数值模拟方法的实现过程。为了方便区分本发明方法和常规方法，以下内容中，将常规隐式差分方法记做isd，将本发明方法记做itsd-1。
图2a至2c为本发明提出的方法与传统方法频散误差对比分析，图中纵坐标为数值模拟速度与真实速度比值，反映了方法的精度。当纵坐标值为1.0时，此时方法无频散。纵坐标偏离1.0的位置越远，则频散越大。图中，速度参数为1500m/s，空间采样步长为15m，算子长度m＝6。图中显示的是传播方向为45度时采用不同时间步长时频散曲线。对比发现，相比传统方法，本发明中提出的方法不同程度的改进了传统方法的精度，压制了数值频散。随着参数n的增大，本发明中提出的方法精度提高，频散变小。
图3为本发明提出的方法与传统方法稳定性对比分析图，图中纵坐标为稳定性因子，当库朗数大于该稳定性因子时，差分方法不稳定。该值越大，方法越稳定。对比发现，传统方法稳定性最差，稳定性因子最小。本发明提出的方法稳定性更好，且随着n的增大，稳定性越来越好。
图4(a)至4(c)为一均匀模型(1200m,2250m)处模拟波场震动图与参考解对比图。模型大小为4500m×4500m，空间步长为15m，时间步长为2ms，采用的子波为主频是18hz的雷克子波，在模型中心激发。参考解为由cagniard-dehoop(dehoop,1960)方法计算得到。本发明中方法采用m＝6。对比发现，传统方法频散误差明显，与参考解吻合差。本发明提出的方法则显著提高了精度，n＝3时可以得到比n＝2更高的精度。图4(d)为各方法与参考解之间的统计均方根误差，明显本发明提出的方法误差要小。
图5为地震勘探数值模拟中常用到的seg/eage盐丘模型。模型大小为5000m×15000m，空间步长为25m，采用的子波为主频是15hz的雷克子波，在地表中心激发，时间步长为2.0ms。本发明中方法采用m＝8和n＝2。图6(a)至6(c)为各方法在2.4s时刻时传播波场快照，参考解为常规显式时空域方法采用很小的时间步长和很大的算子长度计算得到。同参考解对比发现，本发明中方法明显具有更高的精度，白色箭头所指区域频散相比传统方法明显小。图7(a)和(b)为各方法在(500m,4500m)处地震记录震动图相比参考解对比图。图中显示的误差值为统计的均方根误差，明显本发明提出的方法均方根误差更小。
为降低频散误差，传统隐式差分方法采用1ms时间步长时计算的均方根误差为1.1824e-4，消耗的cpu时间为403.58s。本发明提出的方法(itsd-1)采用2ms时间步长时计算的均方根误差为1.2390e-4，消耗的cpu时间为295.27s。可以发现，本发明提出的方法在基本保持了精度的前提下，显著降低了计算时间，效率得到大大提高。
本方法能够显著提高常规模拟方法的精度和效率，能够为计算量巨大的逆时偏移方法和全波形反演方法高效率地提供更为准确的波场信息，具有重要的实际用途。本发明的具体实施方式中未涉及的说明属于本领域公知的技术，可参考公知技术加以实施。
以上具体实施方式及实施例是对本发明提出的一种基于声波方程的线性优化隐式时空域有限差分数值模拟方法技术思想的具体支持，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在本技术方案基础上所做的任何等同变化或等效的改动，均仍属于本发明技术方案保护的范围。