[发明专利]一种基于声波方程的线性优化隐式时空域有限差分数值模拟方法

权利要求书

1.一种基于声波方程的线性优化隐式时空域有限差分数值模拟方法，其特征在于，包括以下几个步骤：

(1)读取参数，所述参数包括正演模拟所需要的速度模型参数文件、有限差分算子长度、子波函数及子波主频、正演所采用的时间与空间步长及地震记录时长；

(2)基于菱形差分算子和二阶时间中心差分，得到时间导数的时间高阶离散格式；

(3)求解二阶空间导数隐式差分格式，在显式高阶差分的基础上，分母中引入二阶中心差分格式，构造空间隐式差分格式求解二阶空间导数，用于减小有限差分算子长度；进一步，基于最小二乘法，以L2范数目标函数极小化空间域频散关系，求解优化的隐式差分系数；

(4)基于所述时间高阶离散格式和空间隐式离散格式，带入声波方程，构建同时具有时间高阶和空间隐式的差分递推格式，得到时空域频散关系；

(5)固定空间差分系数，采用优化策略极小化整个递推格式频散关系，建立关于时间差分格式中差分系数的L2范数目标函数，利用线性优化算法对时间差分系数进行求解；

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

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

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

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

其中，h为空间采样步长，v为速度，a_m,n为有限差分系数；N为菱形算子中的差分算子长度；n和m为菱形算子中的坐标位置序列；

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

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

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

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

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

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

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

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

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

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

其中，ω为角频率，k_x和k_z分别为沿x和z方向的波数，r为库朗数，r＝vτ/h；

步骤(5)中，首先，差分系数a′_m和b由公式(5)计算得到，在计算差分系数a_m,0和a_m,n时保持固定，此时优化差分系数a_m,0和a_m,n是线性优化问题，直接线性求解，定义如下目标函数：

其中，d＝[d₁,L,d_N+(N-1)N/2]^T＝[a_1,0,L,a_N,0,a_1,1,L,a_1,N-1,L,a_N-1,1]^T，表示差分系数；θ为传播角度，q为介于0和1间的实数，波数k_x和k_z与波数k之间满足如下关系：

k_x＝kcos(θ),k_z＝ksin(θ), (13)ψ(β,θ,r)及定义如下：

目标函数是关于差分系数a_m,0和a_m,n的线性函数，差分系数a_m,0和a_m,n通过求解如下公式得到，

2.根据权利要求1所述的基于声波方程的线性优化隐式时空域有限差分数值模拟方法，其特征在于，步骤(6)中，整个计算区域被划分为内部区域、边界区域和过渡区域，具体的方法如下：

(6-1)利用公式(2)和公式(3)离散时间和空间导数，采用公式(5)和公式(16)计算差分系数，进行双程波方程求解，求解声波方程，计算所有区域处的双程波波场；

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

针对二阶声波方程，以左边界和左上角为例，左边界单程波方程为：

其中，x和z代表坐标；

左上角单程波方程为：

其他边界处及角点处单程波可以类似得到，并进行对应求解；采用上述吸收边界条件策略计算各位置处任意时刻波场值；

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

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

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

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