[发明专利]一种频率域地震正演模拟的高精度差分数值法

权利要求书

1.一种频率域地震正演模拟的高精度差分数值法，其特征在于，包括以下步骤：

S1、对二维标量频率域波动方程高阶项离散，首先在25点模板的水平或垂直方向分别进行单倍和双倍步长的二阶离散，然后在垂直或水平方向利用对称性和泰勒展开式消除二阶项，得到两种四阶精度的逼近；

S2、把这两种四阶精度的逼近与经典的四阶中心差分格式进行加权组合；其中，权系数的和为1；

S3、对于低阶项的离散，先利用25点模板四个角处的16个网格点进行离散，得到四阶精度逼近，接着再与水平和竖直方向的四阶中心离散进行加权组合，权系数的和为1；

S4、对得到的差分方案进行频散分析，并基于频散极小化原理得到各个权系数的值。

2.根据权利要求1所述的一种频率域地震正演模拟的高精度差分数值法，其特征在于，已知二维标量频率域波动方程为：

这里，未知量u表示压力场，为波数，其中f和v分别表示波的频率和波的传播速度；

为了在有限的计算区域内进行模拟仿真，对方程(1)进行添加PML边界层，得到

其中C＝e_xe_y，ω为角频率，σ_x定义如下：

其中f₀为震源主频，a₀为一经验常数，L_PML为PML的厚度，l_x为到PML和内部区域间界面的距离；σ_y的定义与σ_x类似；

记u_m,n为函数u在网格点(x_m,y_n)处的离散值，定义符号如下：

其中j∈{-2,-1,1,2},t＝1,2,h为离散步长；定义符号：

继续定义符号Φ_xu为：

其中α₁+α₂+α₃＝1，Φ_xu即为高阶项的离散逼近；

类似地，分别定义符号Φ_yu如下：

Φ_yu即为高阶项的离散逼近；则频率域波动方程(2)的高阶项离散逼近为：

-Φ_xu-Φ_yu (3)

定义符号Ι₁，Ι₂，Ι₃如下，其中j＝{-2,-1,0,1,2}；

I₁(Ck²u)＝(Ck²u)_m,n,

则频率域波动方程(2)的低阶项离散逼近为

I(Ck²u):＝β₁I₁(Ck²u)+β₂Ι₂(Ck²u)+β₃Ι₃(Ck²u), (4)

其中β₁+β₂+β₃＝1，综合(3)式与(4)式，得到频率域波动方程(2)的4阶差分逼近；

-Φ_xu-Φ_yu-I(Ck²u)＝0 (5)

在PML内部区域，25点四阶差分方案(5)简化为

其中，U_m+j,n+l为未知量u在(x_m+j,y_n+l)处的离散值，j,l＝-2,-1,0,1,2，系数T₁至T₆为：

令U(x,y):＝e^{-ik(xcosθ+ysinθ)}，则U(x,y)是频率域波动方程的面波解，其中θ为波传播方向与y轴的夹角；记λ为波长，G为单个波长内的网格点数，即G＝λ/h；又由于λ＝v/f，故有kh＝2π/G；

把离散面波解带入差分格式(6)，并利用欧拉公式e^-ix＝cosx+isinx式进行变换，得到频散方程如下：

4T₁P₂Q₂+4T₂(P₁Q₂+P₂Q₁)+2T₃(P₂+Q₂)+4T₄P₁Q₁+2T₅(P₁+Q₁)+T₄＝0, (7)

其中

把T_j的值带入(7)式，j＝1,2,...,6；并用数值波数k_N替换k，得到

其中

R＝a₁M+a₂N+a₃O,

M＝[(P₂+Q₂)-16(P₁+Q₁)+30],

N＝[2P₂Q₂-4(P₁Q₂+P₂Q₁)-(P₂+Q₂)+4(P₁+Q₁)],

O＝[4(P₁Q₂+P₂Q₁)-4(P₂+Q₂)-32P₁Q₁+16(P₁+Q₁)]

联合kh＝2π/G和(8)式得到

记则求参数α₁,α₂,α₃,β₁,β₂,β₃的值转化为求最优化问题

其中，由对称性得到，而G≥2由Nyqusit采样定理得到；

令并带入α₁＝1-α₂-α₃,β₁＝1-β₂-β₃,得到方程(11)

其中P₁,Q₁,P₂,Q₂,O,M,N是关于θ,G的函数；分别在θ和G的取值范围内采样l＝10和r＝100个点，即

把(12)中的l*r＝1000个点带入方程(11)得到超定线性系统

其中，

利用最小二乘法求解线性系统(13)，即

得到α₂＝0.0266,α₃＝0.1923,β₁＝0.0970,β₂＝0.1902,则α₁＝1-α₂-α₃＝0.7811,

β₁＝1-β₂-β₃＝0.7128

把参数α₁,α₂,α₃,β₁,β₂,β₃的值带入(5)式，得到最终的频率波动方程正演模拟的4阶25点差分方案。