CN106154331A - 正交介质地震波场模拟频散压制方法 - Google Patents

Abstract

本发明提供一种正交介质地震波场模拟频散压制方法，包括：步骤1，利用岩石物理理论和裂隙等效理论将裂隙物性参数转化为正交介质参数，为正交介质正演模拟提供模型；步骤2，根据地震波动力学理论结合裂隙介质刚度矩阵推导正交介质弹性波一阶速度‑应力方程；步骤3，利用有限差分法求解正交介质弹性波方程模拟地震波在裂隙介质中的传播过程；步骤4，利用最小二乘优化方法优化正交介质交错网格差分系数从而压制频散提高模拟精度；步骤5，进行模拟测试。该正交介质地震波场模拟频散压制方法推导更大频率范围内可以获得高精度模拟结果的差分系数，可以用于压制频散，当模型速度低或地震子波主频高时效果更为明显。

Description

正交介质地震波场模拟频散压制方法

技术领域

本发明涉及勘探地球物理领域，特别是涉及到一种正交介质地震波场模拟频散压制方法。

背景技术

地震学研究的对象是地震波及其传播的地球介质，实际地球介质是一种非均匀、非完全弹性、各向异性、多相态的介质。地震学的发展历程正是由简单均匀、完全弹性、各向同性、单相态的波动理论向复杂的真实地球介质的波动理论步步逼近的过程(梁锴,2009)。开展地震各向异性研究对认知地球介质结构、勘探开发复杂油气藏和预报地质灾害等方面均具有理论意义和实用价值，是认知地球介质历史发展的必然(吴国忱，2006)。随着油气勘探开发的深入，勘探的目标逐步由构造油气藏转化为裂缝性油气藏，裂隙是许多油气藏中液体或气体的重要通道，正确的识别裂隙发育位置对于油气藏的勘探和开发均具有重要意义(秦海旭，2015)。裂隙中弹性波传播表现为方位各向异性介质。具有方位各向异性特征的介质包括HTI介质、正交各向异性介质和单斜各向异性介质等。Tsvankin在2001年指出可将发育在横向各向同性介质中的多组平行排列直立裂隙等效为正交各向异性(Tsvankin,2001)，图1给出了正交各向异性介质示意图(Tsvankin,2001)。更多情况下前人研究的是由广泛扩容各向异性(EDA)(Crampin，1987)和旋回性薄互层(PTL)(Postma，1995)等效而成(何燕，2008)。本文研究的是两组正交发育于各向同性介质中的直立裂隙等效成的正交各向异性介质(Bakulin，2000partII)。各向异性介质三大裂隙理论主要包括Hudson裂缝等效理论(Hudson,1981)、Thomsen裂隙等效理论(Thomsen,1995)和线性滑动理论(Schoenberg,1995)。这里利用线性滑动理论将前文提到的裂隙模型等效为正交各向异性来进行岩石物理建模，为地震波数值模拟提供模型。

地震波正演模拟是模拟地震波在地球介质中的传播过程，并研究地震波的传播特性与地球介质参数的关系(孙成禹,2007)，开展地震波的正演模拟研究，对人们正确认识地震波的传播规律，验证所求地球模型的正确性，进行实际地震资料的地质解释与储层预测等，均具有重要的理论和实际意义(吴国忱，2006)。各向异性介质中地震波正演模拟的数值方法主要有：射线追踪法、有限差分法、反射率法、伪谱法和有限元法等，这些数值正演模拟方法各有优缺点。有限差分法是一种应用比较广泛的正演模拟方法，能够较精确地模拟任意非均匀介质中的地震波场，并含有多次散射、转换波与绕射波。有别于规则网格有限差分法，交错网格有限差分法采用一阶速度-应力弹性波方程，其无须对弹性常数进行空间微分(Virieux,1984)，在相同储存空间和计算量的情况下具有更高的模拟精度(董良国，2000)。传统交错网格有限差分空间导数的高阶差分系数一般是通过Taylor级数展开法求取的，而用该方法确定的差分系数来求解空间导数时，一般只是在一个较小的频率范围内才能取得比较高的精度(杨蕾，2014)。刘洋在2013年采用最小二乘法求解二阶空间导数差分系数(Yang Liu,2013)，杨蕾在2014年使用最小二乘法优化一阶空间导数差分系数(杨蕾，2014)。

发明内容

本发明的目的是提供一种可以压制波场模拟中的数值频散现象，明显提高模拟精度的正交介质地震波场模拟频散压制方法。

本发明的目的可通过如下技术措施来实现：正交介质地震波场模拟频散压制方法，该正交介质地震波场模拟频散压制方法包括：步骤1，利用岩石物理理论和裂隙等效理论将裂隙物性参数转化为正交介质参数，为正交介质正演模拟提供模型；步骤2，根据地震波动力学理论结合裂隙介质刚度矩阵推导正交介质弹性波一阶速度-应力方程；步骤3，利用有限差分法求解正交介质弹性波方程模拟地震波在裂隙介质中的传播过程；步骤4，利用最小二乘优化方法优化正交介质交错网格差分系数从而压制频散提高模拟精度；步骤5，进行模拟测试。

本发明的目的还可通过如下技术措施来实现：

在步骤1中，借助线性滑动理论给出了各向同性介质背景下的两组正交直立裂隙的弹性系数矩阵：

$C=\left[\begin{array}{cccccc}{C}_{11}& C_{12}& C_{13}& & & \\ C_{12}& C_{22}& C_{23}& & & \\ C_{13}& C_{23}& C_{33}& & & \\ & & & C_{44}& & \\ & & & & C_{55}& \\ & & & & & C_{66}\end{array}\right]=\left(\begin{array}{cc}{\tilde{C}}_1& 0\\ 0& {\tilde{C}}_2\end{array}\right)---\left(1\right)$

其中C为弹性系数矩阵，C₁₁,C₂₂,C₃₃,C₄₄,C₅₅,C₆₆,C₁₂,C₁₃,C₂₃为弹性系数矩阵元素，为子矩阵，其又可表示为：

${\tilde{C}}_1=\frac{1}{d}\left[\begin{array}{ccc}(\mathrm{\λ}+2\mathrm{\μ})l_1{m}_3& {\mathrm{\λl}}_1{m}_1& {\mathrm{\λl}}_1{m}_2\\ {\mathrm{\λl}}_1{m}_1& (\mathrm{\λ}+2\mathrm{\μ})l_3{m}_1& {\mathrm{\λl}}_2{m}_1\\ {\mathrm{\λl}}_1{m}_2& {\mathrm{\λl}}_2{m}_1& (\mathrm{\λ}+2\mathrm{\μ})(l_1{m}_3-l_4)\end{array}\right]---\left(2\right)$

${\tilde{C}}_2=\left[\begin{array}{ccc}\mathrm{\μ}(1-{\mathrm{\Δ}}_{T2})& 0& 0\\ 0& \mathrm{\μ}(1-{\mathrm{\Δ}}_{T1})& 0\\ 0& 0& \frac{\mathrm{\μ}(1-{\mathrm{\Δ}}_{T1})(1-{\mathrm{\Δ}}_{T2})}{(1-{\mathrm{\Δ}}_{T1}{\mathrm{\Δ}}_{T2})}\end{array}\right]---\left(3\right)$

其中λ和μ分别为裂隙所在背景介质的拉梅参数，Δ_N和Δ_T分别是法向柔度和切向柔度，它们与裂隙充填物有关，变量的下标N1,N2,T1,T2分别表示第一组和第二组裂隙法向和切向，d＝1-r²Δ_N1Δ_N2,r＝1-2g,g＝μ/(λ+2μ)，l₁,l₂,l₃,l₄,m₁,m₂,m₃,d为过渡参数，其可以表示为：

$\begin{array}{c}{l}_1=1-{\mathrm{\Δ}}_{N1},l_2=1-r-{\mathrm{\Δ}}_{N1},l_3=1-r^2{\mathrm{\Δ}}_{N1},l_4=4r^2{g}^2{\mathrm{\Δ}}_{N1}{\mathrm{\Δ}}_{N2}\\ m_1=1-{\mathrm{\Δ}}_{N2},m_2=1-{\mathrm{r\Δ}}_{N2},m_3=1-r^2{\mathrm{\Δ}}_{N2}\end{array}---\left(4\right)$

切向柔量和法向柔量的表达式，其中K为体积模量：

${\mathrm{\Δ}}_{Ni}=\frac{(\mathrm{\λ}+2\mathrm{\μ})K_{Ni}}{1+(\mathrm{\λ}+2\mathrm{\μ})K_{Ni}},{\mathrm{\Δ}}_{Ti}=\frac{{\mathrm{\μK}}_{Ti}}{1+{\mathrm{\μK}}_{Ti}}---\left(5\right)$

当裂隙满足Hudson理论的假设时，可以从Hudson理论中给出Δ_T,Δ_N的表达式：

(1)当裂缝中填充较小体积模量和剪切模量的固体时：

${\mathrm{\Δ}}_N=\frac{4}{3g\[1-g+(k^{\′}+\frac{4{\mathrm{\μ}}^{\′}}{3})/\left(\mathrm{\π}d\mathrm{\μ}\right)\]}e,{\mathrm{\Δ}}_T=\frac{16}{3\[3-2g+\frac{4{\mathrm{\μ}}^{\′}}{\mathrm{\π}d\mathrm{\μ}}\]}e---\left(6\right)$

(2)当裂缝为干裂缝时：

${\mathrm{\Δ}}_N=\frac{4}{3g(1-g)}e,{\mathrm{\Δ}}_T=\frac{16}{3(3-2g)}e---\left(7\right)$

(3)当裂缝中填充无粘滞流体时：

${\mathrm{\Δ}}_N=0,{\mathrm{\Δ}}_T=\frac{16}{3(3-2g)}e---\left(8\right)$

其中，k'和μ'分别为填充介质的体积模量和剪切模量，e和d分别为裂缝密度和裂缝的纵横比。

在步骤2中，根据应力应变关系：

$\left[\begin{array}{c}{\mathrm{\τ}}_{xx}\\ {\mathrm{\τ}}_{yy}\\ {\mathrm{\τ}}_{zz}\\ {\mathrm{\τ}}_{yz}\\ {\mathrm{\τ}}_{xz}\\ {\mathrm{\τ}}_{xy}\end{array}\right]=\left[\begin{array}{cccccc}{c}_{11}& c_{12}& c_{13}& 0& 0& c_{16}\\ c_{12}& c_{22}& c_{23}& 0& 0& c_{26}\\ c_{13}& c_{23}& c_{33}& 0& 0& c_{36}\\ 0& 0& 0& c_{44}& c_{45}& 0\\ 0& 0& 0& c_{45}& c_{55}& 0\\ c_{16}& c_{26}& c_{36}& 0& 0& c_{66}\end{array}\right]\left[\begin{array}{c}{e}_{xx}\\ e_{yy}\\ e_{zz}\\ e_{yz}\\ e_{xz}\\ e_{xy}\end{array}\right]---\left(9\right)$

其中τ_xx,τ_yy,τ_zz,τ_xy,τ_xz,τ_yz为应力张量，e_xx,e_yy,e_zz,e_xy,e_xz,e_yz为应变张量，C₁₁,C₂₂,C₃₃,C₄₄,C₅₅,C₆₆,C₁₂,C₁₃,C₂₃为弹性系数矩阵元素，

应力位移关系：

$\mathrm{\ρ}\frac{{\∂}^2}{\∂t^2}\left[\begin{array}{c}{u}_x\\ u_y\\ u_z\end{array}\right]=\left[\begin{array}{cccccc}\frac{\∂}{\∂x}& 0& 0& 0& \frac{\∂}{\∂z}& \frac{\∂}{\∂y}\\ 0& \frac{\∂}{\∂y}& 0& \frac{\∂}{\∂z}& 0& \frac{\∂}{\∂x}\\ 0& 0& \frac{\∂}{\∂z}& \frac{\∂}{\∂y}& \frac{\∂}{\∂x}& 0\end{array}\right]\left[\begin{array}{c}{\mathrm{\τ}}_{xx}\\ {\mathrm{\τ}}_{yy}\\ {\mathrm{\τ}}_{zz}\\ {\mathrm{\τ}}_{yz}\\ {\mathrm{\τ}}_{xz}\\ {\mathrm{\τ}}_{xy}\end{array}\right]+\mathrm{\ρ}\left[\begin{array}{c}{F}_x\\ F_y\\ F_z\end{array}\right]---\left(10\right)$

其中ρ为密度，u_x,u_y,u_z为位移变量，F_x,F_y,F_z为震源变量，

位移应变关系：

$\left[\begin{array}{c}{e}_{xx}\\ e_{yy}\\ e_{zz}\\ e_{yz}\\ e_{xz}\\ e_{xy}\end{array}\right]=\left[\begin{array}{ccc}\frac{\∂}{\∂x}& 0& 0\\ 0& \frac{\∂}{\∂y}& 0\\ 0& 0& \frac{\∂}{\∂z}\\ 0& \frac{\∂}{\∂z}& \frac{\∂}{\∂y}\\ \frac{\∂}{\∂z}& 0& \frac{\∂}{\∂x}\\ \frac{\∂}{\∂y}& \frac{\∂}{\∂x}& 0\end{array}\right]\left[\begin{array}{c}{u}_x\\ u_y\\ u_z\end{array}\right]---\left(11\right)$

即可推导三维正交介质弹性波一阶速度-应力方程：

$\left\{\begin{array}{c}\mathrm{\ρ}\frac{\∂v_x}{\∂t}=\frac{\∂{\mathrm{\τ}}_{xx}}{\∂x}+\frac{\∂{\mathrm{\τ}}_{xz}}{\∂z}+\frac{\∂{\mathrm{\τ}}_{xy}}{\∂y}+{\mathrm{\ρF}}_x\\ \mathrm{\ρ}\frac{\∂v_y}{\∂t}=\frac{\∂{\mathrm{\τ}}_{yy}}{\∂y}+\frac{\∂{\mathrm{\τ}}_{yz}}{\∂z}+\frac{\∂{\mathrm{\τ}}_{xy}}{\∂x}+{\mathrm{\ρF}}_y\\ \mathrm{\ρ}\frac{\∂v_z}{\∂t}=\frac{\∂{\mathrm{\τ}}_{zz}}{\∂z}+\frac{\∂{\mathrm{\τ}}_{yz}}{\∂y}+\frac{\∂{\mathrm{\τ}}_{xz}}{\∂x}+{\mathrm{\ρF}}_z\end{array}\right.$

$\left\{\begin{array}{c}\frac{\∂{\mathrm{\τ}}_{xx}}{\∂t}=c_{11}\frac{\∂v_x}{\∂x}+c_{12}\frac{\∂v_y}{\∂y}+c_{13}\frac{\∂v_z}{\∂z}+c_{16}(\frac{\∂v_x}{\∂y}+\frac{\∂v_y}{\∂x})\\ \frac{\∂{\mathrm{\τ}}_{yy}}{\∂t}=c_{12}\frac{\∂v_x}{\∂x}+c_{22}\frac{\∂v_y}{\∂y}+c_{23}\frac{\∂v_z}{\∂z}+c_{26}(\frac{\∂v_x}{\∂y}+\frac{\∂v_y}{\∂x})\\ \frac{\∂{\mathrm{\τ}}_{zz}}{\∂t}=c_{13}\frac{\∂v_x}{\∂x}+c_{23}\frac{\∂v_y}{\∂y}+c_{33}\frac{\∂v_z}{\∂z}+c_{36}(\frac{\∂v_x}{\∂y}+\frac{\∂v_y}{\∂x})\\ \frac{\∂{\mathrm{\τ}}_{yz}}{\∂t}=c_{44}\left(\frac{\∂v_y}{\∂z}+\frac{\∂v_z}{\∂y}\right)+c_{45}(\frac{\∂v_x}{\∂z}+\frac{\∂v_z}{\∂x})\\ \frac{\∂{\mathrm{\τ}}_{xz}}{\∂t}=c_{45}\left(\frac{\∂v_y}{\∂z}+\frac{\∂v_z}{\∂y}\right)+c_{55}(\frac{\∂v_x}{\∂z}+\frac{\∂v_z}{\∂x})\\ \frac{\∂{\mathrm{\τ}}_{xy}}{\∂t}=c_{16}\frac{\∂v_x}{\∂x}+c_{26}\frac{\∂v_y}{\∂y}+c_{36}\frac{\∂v_z}{\∂z}+c_{66}(\frac{\∂v_x}{\∂y}+\frac{\∂v_y}{\∂x})\end{array}\right.$

其中，v_x,v_y,v_z为速度变量，t表示时间变量。

在步骤3中，利用泰勒展开得到时间2阶和空间2N阶差分形式：

时间2阶差分：

$\frac{\∂u\left(t\right)}{\∂t}=\frac{u(t+\mathrm{\Δ}t)-u(t-\mathrm{\Δ}t)}{2\mathrm{\Δ}t}---\left(12\right)$

空间2N阶差分：

$\frac{\∂u}{\∂x}=\frac{1}{\mathrm{\Δ}x}\underset{n=1}{\overset{N}{\Σ}}{c}_n\[u(x+(2n-1\left)\frac{\mathrm{\Δ}x}{2}\right)-u(x-(2n-1\left)\frac{\mathrm{\Δ}x}{2}\right)\]+O\left({\mathrm{\Δx}}^{2N}\right)---\left(13\right)$

其中O(Δx^2N)为泰勒展开误差项，N阶差分系数c_m可以通过以下方程求取：

$a\left[\begin{array}{ccccc}{1}^1& 3^1& 5^1& \mathrm{...}& {(2N-1)}^1\\ 1^3& 3^3& 5^3& \mathrm{...}& {(2N-1)}^3\\ 1^5& 3^5& 5^5& \mathrm{...}& {(2N-1)}^5\\ \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}\\ 1^{(2N-1)}& 3^{(2N-1)}& 5^{(2N-1)}& \mathrm{...}& {(2N-1)}^{(2N-1)}\end{array}\right]\left[\begin{array}{c}{c}_1\\ c_2\\ c_3\\ .\\ .\\ .\\ c_N\end{array}\right]=\left[\begin{array}{c}1\\ 0\\ 0\\ .\\ .\\ .\\ 0\end{array}\right]---\left(14\right)$

在步骤4中，将最小二乘优化差分系数的方法应用于三维正交介质模拟中，一阶空间导数Taylor展开：

$\frac{\∂u}{\∂x}\≈\frac{1}{h}\underset{m=1}{\overset{M}{\Σ}}{c}_m\[u(x+mh-0.5h)-u(x-mh+0.5h)\]---\left(15\right)$

其中h为空间采样间隔，c_m为空间差分系数，M为差分阶数，

将平面波解u(x)＝u₀e^ikx带入，其中_k为波数，结合欧拉公式

e^±ix＝cosx±isinx，化简可得：

$\frac{kh}{2}\≈\underset{m=1}{\overset{M}{\Σ}}{a}_m\sin \[kh(2m-1)/2\]---\left(16\right)$

定义过渡函数β和为：

可得：

定义误差函数：

其中E表示误差，a_M为空间差分系数，根据最小二乘基本原理求解：

$\frac{\∂E}{\∂c_m}=0---\left(20\right)$

写成矩阵形式即：

其中c_m为空间差分系数，M为差分阶数，算法和分别为：

本发明中的正交介质地震波场模拟频散压制方法，针对正交介质波场模拟过程中数值频散现象，和地震源子波主频较高或模拟模型速度较低时频散现象尤为明显的情况，借助线性滑动理论，将裂隙密度、纵横比等不同裂隙物性参数等效为正交介质刚度矩阵，提出两组正交直立裂隙介质建模方法。利用交错网格高阶有限差分法实现三维正交介质弹性波波场模拟，并提出正交介质中基于最小二乘法优化的高精度交错网格高阶有限差分法频散压制方法。模拟结果表明介质各向异性强度与裂隙物性直接相关，相比普通交错网格，最小二乘优化方法可以压制波场模拟中的数值频散现象，明显提高模拟精度。

附图说明

图1为正交介质模型的示意图；

图2为本发明的正交介质地震波场模拟频散压制方法的一具体实施例的流程图；

图3为本发明的一具体实施例中正交介质变量定义网格的示意图；

图4为本发明的一具体实施例中优化前后不同差分阶数数值频散随频率变化图；

图5为本发明的一具体实施例中交错网格空间6阶差分均匀介质模型250MS正演波场优化效果的示意图；

图6为本发明的一具体实施例中交错网格空间6阶差分均匀介质模型500MS正演波场优化效果的示意图；

图7为本发明的一具体实施例中交错网格空间8阶差分均匀介质模型500MS正演波场优化效果的示意图。

具体实施方式

为使本发明的上述和其他目的、特征和优点能更明显易懂，下文特举出较佳实施例，并配合附图所示，作详细说明如下。

如图2所示，图2为本发明的正交介质地震波场模拟频散压制方法的流程图。

步骤101：利用岩石物理理论和裂隙等效理论将裂隙物性参数转化为正交介质参数，为正交介质正演模拟提供模型；

该方法首先需要利用岩石物理理论和裂隙介质等效理论建立正交介质模型。Bakulin借助线性滑动理论给出了各向同性介质背景下的两组正交直立裂隙的弹性系数矩阵：

$C=\left[\begin{array}{cccccc}{C}_{11}& C_{12}& C_{13}& & & \\ C_{12}& C_{22}& C_{23}& & & \\ C_{13}& C_{23}& C_{33}& & & \\ & & & C_{44}& & \\ & & & & C_{55}& \\ & & & & & C_{66}\end{array}\right]=\left(\begin{array}{cc}{\tilde{C}}_1& 0\\ 0& {\tilde{C}}_2\end{array}\right)---\left(1\right)$

其中C为弹性系数矩阵，C₁₁,C₂₂,C₃₃,C₄₄,C₅₅,C₆₆,C₁₂,C₁₃,C₂₃为弹性系数矩阵元素，为子矩阵，其又可表示为：

${\tilde{C}}_1=\frac{1}{d}\left[\begin{array}{ccc}(\mathrm{\λ}+2\mathrm{\μ})l_1{m}_3& {\mathrm{\λl}}_1{m}_1& {\mathrm{\λl}}_1{m}_2\\ {\mathrm{\λl}}_1{m}_1& (\mathrm{\λ}+2\mathrm{\μ})l_3{m}_1& {\mathrm{\λl}}_2{m}_1\\ {\mathrm{\λl}}_1{m}_2& {\mathrm{\λl}}_2{m}_1& (\mathrm{\λ}+2\mathrm{\μ})(l_1{m}_3-l_4)\end{array}\right]---\left(2\right)$

${\tilde{C}}_2=\left[\begin{array}{ccc}\mathrm{\μ}(1-{\mathrm{\Δ}}_{T2})& 0& 0\\ 0& \mathrm{\μ}(1-{\mathrm{\Δ}}_{T1})& 0\\ 0& 0& \frac{\mathrm{\μ}(1-{\mathrm{\Δ}}_{T1})(1-{\mathrm{\Δ}}_{T2})}{(1-{\mathrm{\Δ}}_{T1}{\mathrm{\Δ}}_{T2})}\end{array}\right]---\left(3\right)$

其中λ和μ分别为裂隙所在背景介质的拉梅参数，Δ_N和Δ_T分别是法向柔度和切向柔度，它们与裂隙充填物有关，变量的下标N1,N2,T1,T2分别表示第一组和第二组裂隙法向和切向，d＝1-r²Δ_N1Δ_N2,r＝1-2g,g＝μ/(λ+2μ)，l₁,l₂,l₃,l₄,m₁,m₂,m₃,d为过渡参数，其可以表示为：

$\begin{array}{c}{l}_1=1-{\mathrm{\Δ}}_{N1},l_2=1-r-{\mathrm{\Δ}}_{N1},l_3=1-r^2{\mathrm{\Δ}}_{N1},l_4=4r^2{g}^2{\mathrm{\Δ}}_{N1}{\mathrm{\Δ}}_{N2}\\ m_1=1-{\mathrm{\Δ}}_{N2},m_2=1-{\mathrm{r\Δ}}_{N2},m_3=1-r^2{\mathrm{\Δ}}_{N2}\end{array}---\left(4\right)$

切向柔量和法向柔量的表达式，其中K为体积模量：：

${\mathrm{\Δ}}_{Ni}=\frac{(\mathrm{\λ}+2\mathrm{\μ})K_{Ni}}{1+(\mathrm{\λ}+2\mathrm{\μ})K_{Ni}},{\mathrm{\Δ}}_{Ti}=\frac{{\mathrm{\μK}}_{Ti}}{1+{\mathrm{\μK}}_{Ti}}---\left(5\right)$

当裂隙满足Hudson理论的假设时，可以从Hudson理论中给出Δ_T,Δ_N的表达式：

(1)当裂缝中填充较小体积模量和剪切模量的固体时：

${\mathrm{\Δ}}_N=\frac{4}{3g\[1-g+(k^{\′}+\frac{4{\mathrm{\μ}}^{\′}}{3})/\left(\mathrm{\π}d\mathrm{\μ}\right)\]}e,{\mathrm{\Δ}}_T=\frac{16}{3\[3-2g+\frac{4{\mathrm{\μ}}^{\′}}{\mathrm{\π}d\mathrm{\μ}}\]}e---\left(6\right)$

(2)当裂缝为干裂缝时：

${\mathrm{\Δ}}_N=\frac{4}{3g(1-g)}e,{\mathrm{\Δ}}_T=\frac{16}{3(3-2g)}e---\left(7\right)$

(3)当裂缝中填充无粘滞流体时：

${\mathrm{\Δ}}_N=0,{\mathrm{\Δ}}_T=\frac{16}{3(3-2g)}e---\left(8\right)$

其中，k'和μ'分别为填充介质的体积模量和剪切模量，e和d分别为裂缝密度和裂缝的纵横比。

步骤102：根据地震波动力学理论结合裂隙介质刚度矩阵推导正交介质弹性波一阶速度-应力方程；

利用地震波动力学理论结合裂隙各向异性理论推导正交介质介质弹性波一阶速度-应力方程。

$\left[\begin{array}{c}{\mathrm{\τ}}_{xx}\\ {\mathrm{\τ}}_{yy}\\ {\mathrm{\τ}}_{zz}\\ {\mathrm{\τ}}_{yz}\\ {\mathrm{\τ}}_{xz}\\ {\mathrm{\τ}}_{xy}\end{array}\right]=\left[\begin{array}{cccccc}{c}_{11}& c_{12}& c_{13}& 0& 0& c_{16}\\ c_{12}& c_{22}& c_{23}& 0& 0& c_{26}\\ c_{13}& c_{23}& c_{33}& 0& 0& c_{36}\\ 0& 0& 0& c_{44}& c_{45}& 0\\ 0& 0& 0& c_{45}& c_{55}& 0\\ c_{16}& c_{26}& c_{36}& 0& 0& c_{66}\end{array}\right]\left[\begin{array}{c}{e}_{xx}\\ e_{yy}\\ e_{zz}\\ e_{yz}\\ e_{xz}\\ e_{xy}\end{array}\right]---\left(9\right)$

其中τ_xx,τ_yy,τ_zz,τ_xy,τ_xz,τ_yz为应力张量，e_xx,e_yy,e_zz,e_xy,e_xz,e_yz为应变张量。

应力位移关系：

$\mathrm{\ρ}\frac{{\∂}^2}{\∂t^2}\left[\begin{array}{c}{u}_x\\ u_y\\ u_z\end{array}\right]=\left[\begin{array}{cccccc}\frac{\∂}{\∂x}& 0& 0& 0& \frac{\∂}{\∂z}& \frac{\∂}{\∂y}\\ 0& \frac{\∂}{\∂y}& 0& \frac{\∂}{\∂z}& 0& \frac{\∂}{\∂x}\\ 0& 0& \frac{\∂}{\∂z}& \frac{\∂}{\∂y}& \frac{\∂}{\∂x}& 0\end{array}\right]\left[\begin{array}{c}{\mathrm{\τ}}_{xx}\\ {\mathrm{\τ}}_{yy}\\ {\mathrm{\τ}}_{zz}\\ {\mathrm{\τ}}_{yz}\\ {\mathrm{\τ}}_{xz}\\ {\mathrm{\τ}}_{xy}\end{array}\right]+\mathrm{\ρ}\left[\begin{array}{c}{F}_x\\ F_y\\ F_z\end{array}\right]---\left(10\right)$

其中ρ为密度，u_x,u_y,u_z为位移变量，F_x,F_y,F_z为震源变量。

位移应变关系：

$\left[\begin{array}{c}{e}_{xx}\\ e_{yy}\\ e_{zz}\\ e_{yz}\\ e_{xz}\\ e_{xy}\end{array}\right]=\left[\begin{array}{ccc}\frac{\∂}{\∂x}& 0& 0\\ 0& \frac{\∂}{\∂y}& 0\\ 0& 0& \frac{\∂}{\∂z}\\ 0& \frac{\∂}{\∂z}& \frac{\∂}{\∂y}\\ \frac{\∂}{\∂z}& 0& \frac{\∂}{\∂x}\\ \frac{\∂}{\∂y}& \frac{\∂}{\∂x}& 0\end{array}\right]\left[\begin{array}{c}{u}_x\\ u_y\\ u_z\end{array}\right]---\left(11\right)$

即可推导三维正交介质弹性波一阶速度-应力方程：

$\left\{\begin{array}{c}\mathrm{\ρ}\frac{\∂v_x}{\∂t}=\frac{\∂{\mathrm{\τ}}_{xx}}{\∂x}+\frac{\∂{\mathrm{\τ}}_{xz}}{\∂z}+\frac{\∂{\mathrm{\τ}}_{xy}}{\∂y}+{\mathrm{\ρF}}_x\\ \mathrm{\ρ}\frac{\∂v_y}{\∂t}=\frac{\∂{\mathrm{\τ}}_{yy}}{\∂y}+\frac{\∂{\mathrm{\τ}}_{yz}}{\∂z}+\frac{\∂{\mathrm{\τ}}_{xy}}{\∂x}+{\mathrm{\ρF}}_y\\ \mathrm{\ρ}\frac{\∂v_z}{\∂t}=\frac{\∂{\mathrm{\τ}}_{zz}}{\∂z}+\frac{\∂{\mathrm{\τ}}_{yz}}{\∂y}+\frac{\∂{\mathrm{\τ}}_{xz}}{\∂x}+{\mathrm{\ρF}}_z\end{array}\right.$

$\left\{\begin{array}{c}\frac{\∂{\mathrm{\τ}}_{xx}}{\∂t}=c_{11}\frac{\∂v_x}{\∂x}+c_{12}\frac{\∂v_y}{\∂y}+c_{13}\frac{\∂v_z}{\∂z}+c_{16}(\frac{\∂v_x}{\∂y}+\frac{\∂v_y}{\∂x})\\ \frac{\∂{\mathrm{\τ}}_{yy}}{\∂t}=c_{12}\frac{\∂v_x}{\∂x}+c_{22}\frac{\∂v_y}{\∂y}+c_{23}\frac{\∂v_z}{\∂z}+c_{26}(\frac{\∂v_x}{\∂y}+\frac{\∂v_y}{\∂x})\\ \frac{\∂{\mathrm{\τ}}_{zz}}{\∂t}=c_{13}\frac{\∂v_x}{\∂x}+c_{23}\frac{\∂v_y}{\∂y}+c_{33}\frac{\∂v_z}{\∂z}+c_{36}(\frac{\∂v_x}{\∂y}+\frac{\∂v_y}{\∂x})\\ \frac{\∂{\mathrm{\τ}}_{yz}}{\∂t}=c_{44}\left(\frac{\∂v_y}{\∂z}+\frac{\∂v_z}{\∂y}\right)+c_{45}(\frac{\∂v_x}{\∂z}+\frac{\∂v_z}{\∂x})\\ \frac{\∂{\mathrm{\τ}}_{xz}}{\∂t}=c_{45}\left(\frac{\∂v_y}{\∂z}+\frac{\∂v_z}{\∂y}\right)+c_{55}(\frac{\∂v_x}{\∂z}+\frac{\∂v_z}{\∂x})\\ \frac{\∂{\mathrm{\τ}}_{xy}}{\∂t}=c_{16}\frac{\∂v_x}{\∂x}+c_{26}\frac{\∂v_y}{\∂y}+c_{36}\frac{\∂v_z}{\∂z}+c_{66}(\frac{\∂v_x}{\∂y}+\frac{\∂v_y}{\∂x})\end{array}\right.$

其中，v_x,v_y,v_z为速度变量，t表示时间变量。

步骤103：利用有限差分法求解正交介质弹性波方程模拟地震波在裂隙介质中的传播过程；

通常情况多利用泰勒展开可以得到时间2阶和空间2N阶差分形式：

时间2阶差分：

$\frac{\∂u\left(t\right)}{\∂t}=\frac{u(t+\mathrm{\Δ}t)-u(t-\mathrm{\Δ}t)}{2\mathrm{\Δ}t}---\left(12\right)$

空间2N阶差分：

$\frac{\∂u}{\∂x}=\frac{1}{\mathrm{\Δ}x}\underset{n=1}{\overset{N}{\Σ}}{c}_n\[u(x+(2n-1\left)\frac{\mathrm{\Δ}x}{2}\right)-u(x-(2n-1\left)\frac{\mathrm{\Δ}x}{2}\right)\]+O\left({\mathrm{\Δx}}^{2N}\right)---\left(13\right)$

其中O(Δx^2N)为泰勒展开误差项，方程不同变量在三维有限差分网格中的位置如图3和表1所示。其中差分系数可以通过以下方程求取，如表2所示：

$\left[\begin{array}{ccccc}{1}^1& 3^1& 5^1& \mathrm{...}& {(2N-1)}^1\\ 1^3& 3^3& 5^3& \mathrm{...}& {(2N-1)}^3\\ 1^5& 3^5& 5^5& \mathrm{...}& {(2N-1)}^5\\ \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}\\ 1^{(2N-1)}& 3^{(2N-1)}& 5^{(2N-1)}& \mathrm{...}& {(2N-1)}^{(2N-1)}\end{array}\right]\left[\begin{array}{c}{c}_1\\ c_2\\ c_3\\ .\\ .\\ .\\ c_N\end{array}\right]=\left[\begin{array}{c}1\\ 0\\ 0\\ .\\ .\\ .\\ 0\end{array}\right]---\left(14\right)$

表1正交介质弹性参数变量的空间位置表

表2普通交错网格差分系数

步骤104：利用最小二乘优化方法优化正交介质交错网格差分系数从而压制频散提高模拟精度；

将最小二乘优化差分系数的方法应用于三维正交介质模拟中。一阶空间导数Taylor展开：

$\frac{\∂u}{\∂x}\≈\frac{1}{h}\underset{m=1}{\overset{M}{\Σ}}{c}_m\[u(x+mh-0.5h)-u(x-mh+0.5h)\]---\left(15\right)$

其中h为空间采样间隔，c_m为空间差分系数，M为差分阶数，

将平面波解u(x)＝u₀e^ikx带入，其中k为波数，结合欧拉公式e^±ix＝cosx±isinx，化简可得：

$\frac{kh}{2}\≈\underset{m=1}{\overset{M}{\Σ}}{a}_m\sin \[kh(2m-1)/2\]---\left(16\right)$

定义过渡函数β和为：

可得：

定义误差函数：

其中E表示误差，a_M为空间差分系数，根据最小二乘基本原理求解：

$\frac{\∂E}{\∂c_m}=0---\left(20\right)$

写成矩阵形式即：

其中c_m为空间差分系数，M为差分阶数，算法和分别为：

若给定d＝0.75可求解新的差分系数如表3所示：

表3优化交错网格差分系数

步骤105：进行模拟测试，模型测试说明最小二乘数值频散压制技术对频散的压制效果和模拟精度的提高效果。

图4展示了不同子波主频、不同差分阶数情况下优化差分系数前后理论模拟精度对比，由图我们可以看出不同情况下模拟差分精度均有所提高。图5、6、7展示了不同差分阶数和不同旅行时差分系数优化前后正交介质正演波场情况，从图中可以看出最小二乘优化差分方法可以有效提高差分精度压制数值频散。

本发明借助线性滑动理论，将裂隙密度、纵横比等不同裂隙物性参数等效为正交介质刚度矩阵，提出两组正交直立裂隙介质建模方法。利用交错网格高阶有限差分法实现三维正交介质弹性波波场模拟，并提出正交介质中基于最小二乘法优化的高精度交错网格高阶有限差分法频散压制方法,使用最小二乘方法优化正交介质一阶速度-应力方程空间导数的差分系数，推导更大频率范围内可以获得高精度模拟结果的差分系数。模拟结果表明介质各向异性强度与裂隙物性直接相关，相比普通交错网格，最小二乘优化方法可以压制波场模拟中的数值频散现象，明显提高模拟精度，子波主频较高和模型波速较低时改善效果更为突出。

Claims (5)

1.正交介质地震波场模拟频散压制方法，其特征在于，该正交介质地震波场模拟频散压制方法包括：

步骤1，利用岩石物理理论和裂隙等效理论将裂隙物性参数转化为正交介质参数，为正交介质正演模拟提供模型；

步骤2，根据地震波动力学理论结合裂隙介质刚度矩阵推导正交介质弹性波一阶速度-应力方程；

步骤3，利用有限差分法求解正交介质弹性波方程模拟地震波在裂隙介质中的传播过程；

步骤4，利用最小二乘优化方法优化正交介质交错网格差分系数从而压制频散提高模拟精度；

步骤5，进行模拟测试。

2.根据权利要求1所述的正交介质地震波场模拟频散压制方法，其特征在于，在步骤1中，借助线性滑动理论给出了各向同性介质背景下的两组正交直立裂隙的弹性系数矩阵：

$C=\left[\begin{array}{cccccc}{C}_{11}& C_{12}& C_{13}& & & \\ C_{12}& C_{22}& C_{23}& & & \\ C_{13}& C_{23}& C_{33}& & & \\ & & & C_{44}& & \\ & & & & C_{55}& \\ & & & & & C_{66}\end{array}\right]=\left(\begin{array}{cc}{\tilde{C}}_1& 0\\ 0& {\tilde{C}}_2\end{array}\right)---\left(1\right)$

其中C为弹性系数矩阵，C₁₁,C₂₂,C₃₃,C₄₄,C₅₅,C₆₆,C₁₂,C₁₃,C₂₃为弹性系数矩阵元素，为子矩阵，其又可表示为：

${\tilde{C}}_1=\frac{1}{d}\left[\begin{array}{ccc}(\mathrm{\λ}+2\mathrm{\μ})l_1{m}_3& {\mathrm{\λl}}_1{m}_1& {\mathrm{\λl}}_1{m}_2\\ {\mathrm{\λl}}_1{m}_1& (\mathrm{\λ}+2\mathrm{\μ})l_3{m}_1& {\mathrm{\λl}}_2{m}_1\\ {\mathrm{\λl}}_1{m}_2& {\mathrm{\λl}}_2{m}_1& (\mathrm{\λ}+2\mathrm{\μ})(l_1{m}_3-l_4)\end{array}\right]---\left(2\right)$

${\tilde{C}}_2=\left[\begin{array}{ccc}\mathrm{\μ}(1-{\mathrm{\Δ}}_{T2})& 0& 0\\ 0& \mathrm{\μ}(1-{\mathrm{\Δ}}_{T1})& 0\\ 0& 0& \frac{\mathrm{\μ}(1-{\mathrm{\Δ}}_{T1})(1-{\mathrm{\Δ}}_{T2})}{(1-{\mathrm{\Δ}}_{T1}{\mathrm{\Δ}}_{T2})}\end{array}\right]---\left(3\right)$

其中λ和μ分别为裂隙所在背景介质的拉梅参数，Δ_N和Δ_T分别是法向柔度和切向柔度，它们与裂隙充填物有关，变量的下标N1,N2,T1,T2分别表示第一组和第二组裂隙法向和切向，d＝1-r²Δ_N1Δ_N2,r＝1-2g,g＝μ/(λ+2μ)，l₁,l₂,l₃,l₄,m₁,m₂,m₃,d为过渡参数，其可以表示为：

$\begin{array}{c}{l}_1=1-{\mathrm{\Δ}}_{N1},l_2=1-r-{\mathrm{\Δ}}_{N1},l_3=1-r^2{\mathrm{\Δ}}_{N1},l_4=4r^2{g}^2{\mathrm{\Δ}}_{N1}{\mathrm{\Δ}}_{N2}\\ m_1=1-{\mathrm{\Δ}}_{N2},m_2=1-{\mathrm{r\Δ}}_{N2},m_3=1-r^2{\mathrm{\Δ}}_{N2}\end{array}---\left(4\right)$

切向柔量和法向柔量的表达式，其中K为体积模量：

${\mathrm{\Δ}}_{Ni}=\frac{(\mathrm{\λ}+2\mathrm{\μ})K_{Ni}}{1+(\mathrm{\λ}+2\mathrm{\μ})K_{Ni}},{\mathrm{\Δ}}_{Ti}=\frac{{\mathrm{\μK}}_{Ti}}{1+{\mathrm{\μK}}_{Ti}}---\left(5\right)$

当裂隙满足Hudson理论的假设时，可以从Hudson理论中给出Δ_T,Δ_N的表达式：

(1)当裂缝中填充较小体积模量和剪切模量的固体时：

${\mathrm{\Δ}}_N=\frac{4}{3g\[1-g+(k^{\′}+\frac{4{\mathrm{\μ}}^{\′}}{3})/\left(\mathrm{\π}d\mathrm{\μ}\right)\]}e,{\mathrm{\Δ}}_T=\frac{16}{3\[3-2g+\frac{4{\mathrm{\μ}}^{\′}}{\mathrm{\π}d\mathrm{\μ}}\]}e---\left(6\right)$

(2)当裂缝为干裂缝时：

${\mathrm{\Δ}}_N=\frac{4}{3g(1-g)}e,{\mathrm{\Δ}}_T=\frac{16}{3(3-2g)}e---\left(7\right)$

(3)当裂缝中填充无粘滞流体时：

Δ_N＝0，

其中，k'和μ'分别为填充介质的体积模量和剪切模量，e和d分别为裂缝密度和裂缝的纵横比。

3.根据权利要求1所述的正交介质地震波场模拟频散压制方法，其特征在于，在步骤2中，根据应力应变关系：

$\left[\begin{array}{c}{\mathrm{\τ}}_{xx}\\ {\mathrm{\τ}}_{yy}\\ {\mathrm{\τ}}_{zz}\\ {\mathrm{\τ}}_{yz}\\ {\mathrm{\τ}}_{xz}\\ {\mathrm{\τ}}_{xy}\end{array}\right]=\left[\begin{array}{cccccc}{c}_{11}& c_{12}& c_{13}& 0& 0& c_{16}\\ c_{12}& c_{22}& c_{23}& 0& 0& c_{26}\\ c_{13}& c_{23}& c_{33}& 0& 0& c_{36}\\ 0& 0& 0& c_{44}& c_{45}& 0\\ 0& 0& 0& c_{45}& c_{55}& 0\\ c_{16}& c_{26}& c_{36}& 0& 0& c_{66}\end{array}\right]\left[\begin{array}{c}{e}_{xx}\\ e_{yy}\\ e_{zz}\\ e_{yz}\\ e_{xz}\\ e_{xy}\end{array}\right]---\left(9\right)$

其中τ_xx,τ_yy,τ_zz,τ_xy,τ_xz,τ_yz为应力张量，e_xx,e_yy,e_zz,e_xy,e_xz,e_yz为应变张量，C₁₁,C₂₂,C₃₃,C₄₄,C₅₅,C₆₆,C₁₂,C₁₃,C₂₃为弹性系数矩阵元素，

应力位移关系：

$\mathrm{\ρ}\frac{{\∂}^2}{\∂t^2}\left[\begin{array}{c}{u}_x\\ u_y\\ u_z\end{array}\right]=\left[\begin{array}{cccccc}\frac{\∂}{\∂x}& 0& 0& 0& \frac{\∂}{\∂z}& \frac{\∂}{\∂y}\\ 0& \frac{\∂}{\∂y}& 0& \frac{\∂}{\∂z}& 0& \frac{\∂}{\∂x}\\ 0& 0& \frac{\∂}{\∂z}& \frac{\∂}{\∂y}& \frac{\∂}{\∂x}& 0\end{array}\right]\left[\begin{array}{c}{\mathrm{\τ}}_{xx}\\ {\mathrm{\τ}}_{yy}\\ {\mathrm{\τ}}_{zz}\\ {\mathrm{\τ}}_{yz}\\ {\mathrm{\τ}}_{xz}\\ {\mathrm{\τ}}_{xy}\end{array}\right]+\mathrm{\ρ}\left[\begin{array}{c}{F}_x\\ F_y\\ F_z\end{array}\right]---\left(10\right)$

其中ρ为密度，u_x,u_y,u_z为位移变量，F_x,F_y,F_z为震源变量，

位移应变关系：

$\left[\begin{array}{c}{e}_{xx}\\ e_{yy}\\ e_{zz}\\ e_{yz}\\ e_{xz}\\ e_{xy}\end{array}\right]=\left[\begin{array}{ccc}\frac{\∂}{\∂x}& 0& 0\\ 0& \frac{\∂}{\∂y}& 0\\ 0& 0& \frac{\∂}{\∂z}\\ 0& \frac{\∂}{\∂z}& \frac{\∂}{\∂y}\\ \frac{\∂}{\∂z}& 0& \frac{\∂}{\∂x}\\ \frac{\∂}{\∂y}& \frac{\∂}{\∂x}& 0\end{array}\right]\left[\begin{array}{c}{u}_x\\ u_y\\ u_z\end{array}\right]---\left(11\right)$

即可推导三维正交介质弹性波一阶速度-应力方程：

$\left\{\begin{array}{c}\mathrm{\ρ}\frac{\∂v_x}{\∂t}=\frac{\∂{\mathrm{\τ}}_{xx}}{\∂x}+\frac{\∂{\mathrm{\τ}}_{xz}}{\∂z}+\frac{\∂{\mathrm{\τ}}_{xy}}{\∂y}+{\mathrm{\ρF}}_x\\ \mathrm{\ρ}\frac{\∂v_y}{\∂t}=\frac{\∂{\mathrm{\τ}}_{yy}}{\∂y}+\frac{\∂{\mathrm{\τ}}_{yz}}{\∂z}+\frac{\∂{\mathrm{\τ}}_{xy}}{\∂x}+{\mathrm{\ρF}}_y\\ \mathrm{\ρ}\frac{\∂v_z}{\∂t}=\frac{\∂{\mathrm{\τ}}_{zz}}{\∂z}+\frac{\∂{\mathrm{\τ}}_{yz}}{\∂y}+\frac{\∂{\mathrm{\τ}}_{xz}}{\∂x}+{\mathrm{\ρF}}_z\end{array}\right.$

$\left\{\begin{array}{c}\frac{\∂{\mathrm{\τ}}_{xx}}{\∂t}=c_{11}\frac{\∂v_x}{\∂x}+c_{12}\frac{\∂v_y}{\∂y}+c_{13}\frac{\∂v_z}{\∂z}+c_{16}(\frac{\∂v_x}{\∂y}+\frac{\∂v_y}{\∂x})\\ \frac{\∂{\mathrm{\τ}}_{yy}}{\∂t}=c_{12}\frac{\∂v_x}{\∂x}+c_{22}\frac{\∂v_y}{\∂y}+c_{23}\frac{\∂v_z}{\∂z}+c_{26}(\frac{\∂v_x}{\∂y}+\frac{\∂v_y}{\∂x})\\ \frac{\∂{\mathrm{\τ}}_{zz}}{\∂t}=c_{13}\frac{\∂v_x}{\∂x}+c_{23}\frac{\∂v_y}{\∂y}+c_{33}\frac{\∂v_z}{\∂z}+c_{36}(\frac{\∂v_x}{\∂y}+\frac{\∂v_y}{\∂x})\\ \frac{\∂{\mathrm{\τ}}_{yz}}{\∂t}=c_{44}\left(\frac{\∂v_y}{\∂z}+\frac{\∂v_z}{\∂y}\right)+c_{45}(\frac{\∂v_x}{\∂z}+\frac{\∂v_z}{\∂x})\\ \frac{\∂{\mathrm{\τ}}_{xz}}{\∂t}=c_{45}\left(\frac{\∂v_y}{\∂z}+\frac{\∂v_z}{\∂y}\right)+c_{55}(\frac{\∂v_x}{\∂z}+\frac{\∂v_z}{\∂x})\\ \frac{\∂{\mathrm{\τ}}_{xy}}{\∂t}=c_{16}\frac{\∂v_x}{\∂x}+c_{26}\frac{\∂v_y}{\∂y}+c_{36}\frac{\∂v_z}{\∂z}+c_{66}(\frac{\∂v_x}{\∂y}+\frac{\∂v_y}{\∂x})\end{array}\right.$

其中，v_x,v_y,v_z为速度变量，t表示时间变量。

4.根据权利要求1所述的正交介质地震波场模拟频散压制方法，其特征在于，在步骤3中，利用泰勒展开得到时间2阶和空间2N阶差分形式：

时间2阶差分：

$\frac{\∂u\left(t\right)}{\∂t}=\frac{u(t+\mathrm{\Δ}t)-u(t-\mathrm{\Δ}t)}{2\mathrm{\Δ}t}---\left(12\right)$

空间2N阶差分：

$\frac{\∂u}{\∂x}=\frac{1}{\mathrm{\Δ}x}\underset{n=1}{\overset{N}{\Σ}}{c}_n\[u(x+(2n-1\left)\frac{\mathrm{\Δ}x}{2}\right)-u(x-(2n-1\left)\frac{\mathrm{\Δ}x}{2}\right)\]+O\left({\mathrm{\Δx}}^{2N}\right)---\left(13\right)$

其中O(Δx^2N)为泰勒展开误差项，N阶差分系数c_m可以通过以下方程求取：

$\left[\begin{array}{ccccc}{1}^1& 3^1& 5^1& \mathrm{...}& {(2N-1)}^1\\ 1^3& 3^3& 5^3& \mathrm{...}& {(2N-1)}^3\\ 1^5& 3^5& 5^5& \mathrm{...}& {(2N-1)}^5\\ \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}\\ 1^{(2N-1)}& 3^{(2N-1)}& 5^{(2N-1)}& \mathrm{...}& {(2N-1)}^{(2N-1)}\end{array}\right]\left[\begin{array}{c}{c}_1\\ c_2\\ c_3\\ .\\ .\\ .\\ c_N\end{array}\right]=\left[\begin{array}{c}1\\ 0\\ 0\\ .\\ .\\ .\\ 0\end{array}\right]---\left(14\right).$

5.根据权利要求1所述的正交介质地震波场模拟频散压制方法，其特征在于，在步骤4中，将最小二乘优化差分系数的方法应用于三维正交介质模拟中，一阶空间导数Taylor展开：

$\frac{\∂u}{\∂x}\≈\frac{1}{h}\underset{m=1}{\overset{M}{\Σ}}{c}_m\[u(x+mh-0.5h)-u(x-mh+0.5h)\]---\left(15\right)$

其中h为空间采样间隔，c_m为空间差分系数，M为差分阶数，

将平面波解u(x)＝u₀e^ikx带入，其中_k为波数，结合欧拉公式

e^±ix＝cos x±i sin x，化简可得：

$\frac{kh}{2}\≈\underset{m=1}{\overset{M}{\Σ}}{a}_m\sin \[kh(2m-1)/2\]---\left(16\right)$

定义过渡函数β和为：

可得：

定义误差函数：

其中E表示误差，a_M为空间差分系数，根据最小二乘基本原理求解：

$\frac{\∂E}{\∂c_m}=0---\left(20\right)$

写成矩阵形式即：

其中c_m为空间差分系数，M为差分阶数，算法和分别为：