0 引言

高效、高精度的地震波场数值模拟是全波形反演和逆时偏移成像的基础^[1-2]。理论和实践证明，相较于有限元法、伪谱法、谱元法、边界元法以及无网格法等数值计算方法，有限差分法具有内存占用小、计算效率高的优点，因而被广泛应用。

提高有限差分格式计算精度的有效方法之一是增加计算的节点数，构造更高阶差分格式，但同时运算量和存储需求也随之增加。紧致有限差分(Compact Difference，CD)格式是一种隐式有限差分格式，与显式格式相比，当使用相同数量的计算节点时，计算精度更高，并且对于边界条件的限制较小。CD格式最早用于求解复杂的流体力学问题，Hirsh^[3]基于Hermite公式提出了一种空间四阶的CD格式，以解决求解抛物方程的巨大运算量问题；Dennis等^[4]将空间四阶CD格式应用于求解流扩散方程，使用较粗的空间网格获得了较高的模拟精度；Lele^[5]提出了具有谱方法精度的对称型CD格式，并给出了求解一阶和二阶导数的三对角和五对角形式。引入迎风机制可有效提高CD格式的稳定性，刘宏等^[6]提出了一种具有三阶精度的迎风型CD格式，并用于求解Naiver-Stokes方程；Fu等^[7]为处理多尺度复杂流体的数值计算问题，构造了五点五阶的迎风型CD格式；Chu等^[8]提出了具有六阶精度的三点组合型CD格式，并将该格式拓展到了交错网格数值模拟^[9]；杨宽德^[10]等发展了基于通量校正输运的CD格式，有效解决了粗网格条件下模拟波场的数值频散和震源噪声问题。

在地震波数值模拟领域，王书强等^[11]采用空间五点六阶的CD格式实现了弹性波的数值模拟；Du等^[12]研究了交错网格紧致差分格式，并模拟了横向各向同性介质中的弹性波场；Liao^[13]研究了时间四阶、空间四阶的对称型CD格式，并给出了三维声波方程的数值模拟结果；Zhang等^[14]结合单程波动方程偏移中的系数优化和分裂思想提出了高阶CD格式的分裂算法，并将其应用于声波方程数值模拟；汪勇等^[15-16]利用优化差分系数的紧致差分格式求解二维声波方程和黏滞声波方程；Liao等^[17]基于Padé近似，采用时间和空间四阶CD格式求解声波方程，获得了高精度的声波波场。

应用CD格式求解波动方程时，涉及到计算和存储大型稀疏对角矩阵方程组。当模型网格数较大时，计算量和存储量会大幅上升。提高数值模拟效率的有效策略之一是应用可变网格思想，针对不同速度区域采用不同网格步长进行求解。传统的变网格方法主要通过改变不同网格间的差分系数或在粗、细网格间的过渡区域应用插值技术实现变网格处理。张慧等^[18]将局部变时间步长和变空间网格步长相结合，发展了时空双变网格算法；Huang等^[19]、李振春等^[20]和曲英铭等^[21]分别将双变网格算法应用于起伏地表正演模拟、逆时偏移成像以及全波形反演等领域。可变网格思想具有良好的应用前景，但传统的变网格算法仍存在以下问题：①当参数设置不合理时，传统方法在粗、细网格的过渡区域容易出现虚假反射^[22]，严重影响模拟和成像精度；②仅能针对单一目标进行网格细化，无法面向多个目标区域；③难以实现自动化网格离散过程，需要人工进行网格剖分，且网格划分过程较复杂。Wang等^[23]提出了一种根据模型速度自适应调节网格步长大小的变网格方法，可有效解决上述问题，该方法无须设置过渡带或改变差分系数，能够高效、准确地模拟地震波场。

为了能够准确、高效地进行数值模拟和偏移成像，本文首先构造了求解二维声波方程的时间二阶中心差分、空间四阶紧致差分的离散格式；然后讨论了该格式的计算精度和稳定性条件；为提高计算效率，引入自适应变网格策略，推导了修正声波方程，并应用于逆时偏移成像；最后，利用层状模型和Mar-mousi模型验证了本文方法的正确性和优越性。

1 紧致有限差分原理

求解二阶导数的CD格式可表示为

$ \begin{aligned} & \beta f_{i-2}^{\prime \prime}+\alpha f_{i-1}^{\prime \prime}+f_i^{\prime \prime}+\alpha f_{i+1}^{\prime \prime}+\beta f_{i+2}^{\prime \prime} \\ &= c \frac{f_{i+3}-2 f_i+f_{i-3}}{9 h^2}+b \frac{f_{i+2}-2 f_i+f_{i-2}}{4 h^2}+ \\ & a \frac{f_{i+1}-2 f_i+f_{i-1}}{h^2} \end{aligned} $

(1)

式中：f_i和f″_i分别表示节点i处变量f的函数值和二阶导数值；h为空间网格间距；α、β、a、b和c为差分系数。利用泰勒展开求解差分系数，可以得到一种求解二阶导数的四阶CD(CD4)格式

$ \alpha f_{i-1}^{\prime \prime}+f_i^{\prime \prime}+\alpha f_{i+1}^{\prime \prime}=a \frac{f_{i+1}-2 f_i+f_{i-1}}{h^2} $

(2)

式中：α=0.1；a=1.2。由式(1)和式(2)可以看出，当采用CD格式近似函数在某点的二阶导数值时，不仅需要函数在该点及相邻节点处的函数值，还需要该点及相邻节点处的二阶导数值，因此，CD格式需使用递推法进行数值求解。此外，采用CD4格式近似二阶导数值时，需要三个计算节点以及一层边界条件，而四阶中心差分(FD4)格式需要五个计算节点以及两层边界条件，六阶中心差分(FD6)格式则需要七个计算节点以及三层边界条件，因此，CD4格式对边界条件的要求较低。

对式(2)进行Fourier变换，并利用欧拉公式化简，可得

$ \tilde{k} h=\sqrt{\frac{2 a(1-\cos w)}{2 \alpha \cos w+1}} $

(3)

式中：w=kh，为真波数k与空间步长h的乘积，其取值范围为(0，π)；$\tilde{k}$表示数值波数。定义$\widetilde{w}=\tilde{k} h$，若利用式(2)近似函数的二阶导数时不存在误差，则有$\widetilde{w}=w$。$\widetilde{w}$与w的差别越小，说明差分格式的精度越高，压制数值频散的能力越强，反之则说明该格式的精度较低。

同理，对于2M阶中心差分格式，有

$ \widetilde{w}=\tilde{k} h=\sqrt{2 \sum\limits_{m=1}^M C_m^{(M)}[1-\cos (m w)]} $

(4)

式中$C_m^{(M)}$为2M阶中心差分格式差分系数。

为便于精度比较和分析，以w为横轴，以$\widetilde{w}^2$为纵轴，绘制了CD4、FD4和FD6格式的误差曲线(图 1)。由图 1可以看出，CD4格式的误差曲线与FD6格式较为接近，并显著优于FD4格式。

2 声波方程的紧致差分离散格式

考虑各向同性介质中二维声波方程

$ \frac{1}{v^2} \frac{\partial^2 u}{\partial t^2}=\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial z^2} $

(5)

式中：u为位移场；v为速度场。

采用时间二阶中心差分、空间CD4格式对式(5)进行离散，有

$ \begin{gathered} u_{i, j}^{n+1}=2 u_{i, j}^n-u_{i, j}^{n-1}+ \\ \frac{v^2(\Delta t)^2}{(\Delta x)^2}\left(u_{*, j}^n\right)^{\prime \prime}+\frac{v^2(\Delta t)^2}{(\Delta z)^2}\left(u_{i, *}^n\right)^{\prime \prime} \end{gathered} $

(6)

式中：$u_{i, j}^n$表示n时刻节点(i，j)处的位移；$\left(u_{*, j}^n\right)^{\prime \prime}$和$\left(u_{i, *}^n\right)^{\prime \prime}$分别为∂²u/∂x²和∂²u/∂z²的离散形式；Δt为时间步长；Δx和Δz分别为x和z方向的空间采样间隔。以z方向为例，对空间二阶偏导数项的求解可表示为矩阵形式

$ \left[\begin{array}{ccccc} 1 & \alpha & & & \\ \alpha & 1 & \alpha & & \\ & \ddots & \ddots & \ddots & \\ & & \alpha & 1 & \alpha \\ & & & \alpha & 1 \end{array}\right]_{J \times J}\left[\begin{array}{c} u_{i, 1}^{\prime \prime} \\ u_{i, 2}^{\prime \prime} \\ \vdots \\ u_{i, J-1}^{\prime \prime} \\ u_{i, J}^{\prime \prime} \end{array}\right]=\frac{1}{(\Delta z)^2}\left[\begin{array}{ccccc} -2 a & a & & & \\ a & -2 a & a & & \\ & \ddots & \ddots & \ddots & \\ & & a & -2 a & a \\ & & & a & -2 a \end{array}\right]_{J \times J}\left[\begin{array}{c} u_{i, 1} \\ u_{i, 2} \\ \vdots \\ u_{i, J-1} \\ u_{i, J} \end{array}\right] $

(7)

式中J为z方向的样点数。上式为典型的三对角矩阵方程组，可简记为

$ \boldsymbol{B}_J \boldsymbol{F}_i=\frac{1}{(\Delta z)^2} \boldsymbol{A}_J \boldsymbol{f}_i $

(8)

同理，对于x方向，有

$ \boldsymbol{B}_I \boldsymbol{F}_j=\frac{1}{(\Delta x)^2} \boldsymbol{A}_I \boldsymbol{f}_j $

(9)

式中I为x方向的样点数。

式(8)和式(9)可直接求逆或采用追赶法求解，本文采用计算速度更快的后者。

应用平面波理论对式(6)进行时空域频散分析，进一步分析其适用条件。令

$ u_{i, j}^n=\exp \left[\mathrm{i}\left(i \Delta x k_x+j \Delta z k_z-k v_{\text {num }} n \Delta t\right)\right] $

(10)

式中：$\mathrm{i}=\sqrt{-1}$；k_x=kcosθ、k_z=ksinθ分别为x方向和z方向的波数，θ为波的传播方向与x轴之间的夹角；v_num为数值传播速度。不失一般性，令Δx=Δz=h，将式(10)代入式(6)，化简可得频散关系方程

$ \begin{array}{l} \cos \left(k v_{\text {num }} \Delta t\right)=1-r^2\left\{\frac{a[1-\cos (w \cos \theta)]}{2 \alpha \cos (w \cos \theta)+1}+\right. \\ \left.\frac{a[1-\cos (w \sin \theta)]}{2 \alpha \cos (w \sin \theta)+1}\right\}=g(w, \theta) \end{array} $

(11)

式中r=vΔt/h为Courant数。定义数值波速v_num与真实速度v的比为

$ \gamma=\frac{v_{\text {num }}}{v}=\frac{k v_{\text {num }} \Delta t}{k v \Delta t}=\frac{k v_{\text {num }} \Delta t}{\frac{v \Delta t}{h} k h}=\frac{k v_{\text {num }} \Delta t}{r w} $

(12)

利用式(11)可求取kv_numΔt，并代入式(12)，有

$ \gamma=\frac{\arccos g(w, \theta)}{r w} $

(13)

理想情况下，若无频散，则γ=1；γ与1的偏离程度越大，说明该格式的频散越严重，反之则说明该格式能够有效压制数值频散。取w (0，π)为横轴，γ为纵轴，由于单位波长内的样点数可表示为2π/w，因此横轴可以看作单位波长内的采样点数由∞逐渐减小到2。图 2为r和θ变化时，CD4、FD4和FD6格式的γ曲线。从图 2可以看出：①随单位波长内采样点数逐渐减少，三种方法的数值频散逐渐加剧，且CD4和FD6格式明显优于FD4格式；②当θ取值为0和π/8时，CD4与FD6格式的速度比曲线较为接近，但后者略优于前者；当θ取值为π/4时，CD4和FD6格式的速度比曲线差异增大；③当r取0.3、θ取π/4时，三种格式的速度比曲线均表现出“上翘”特征，且在w∈(2π/5，4π/5)的部分区间内，CD4格式的γ曲线较FD6格式更接近于1。这表明，在有些情况下，CD4格式优于FD6格式。

稳定性条件是衡量时空域差分格式应用效果的重要问题。对式(6)进行Fourier分析，可得该式的稳定性条件

$ v \Delta t \sqrt{\frac{1}{(\Delta x)^2}+\frac{1}{(\Delta z)^2}} \leqslant \sqrt{\frac{1-2 \alpha}{a}} \approx 0.816 $

(14)

在正演模拟过程中，需要对模型外加人工边界，对传播至边界的波场进行吸收、衰减处理^[24-25]。为提高计算效率，本文采用计算量较小的海绵层吸收边界条件，其衰减因子G为

$ \left\{\begin{array}{l} G(i)=\exp \left\{-[q(i-d)]^2\right\} \\ q=\sqrt{-\frac{\lg 0.9}{d^2}} \end{array}\right. $

(15)

式中：q为衰减系数；d为海绵吸收层厚度；i为计算点到外层边界的距离。

3 自适应变网格基本原理

采用式(6)进行声波方程延拓时，计算量主要源于求解三对角矩阵方程，随着模型网格数的增加，计算量也随之增大。为减小计算量，同时保证计算精度，可引入自适应变网格策略，通过优化网格剖分减少模型网格数，从而提高计算效率。

给定网格点数和网格间距已知的规则网格模型，重新计算理论最优的垂向网格间距

$ \Delta z(z)=\frac{v_{\min }(z)}{l F_{\max }} $

(16)

式中：v_min(z)为模型纵向上每一层的最小速度；F_max表示震源子波的最大频率；l为每个波长内的采样点数，一般取6或7，为避免数值频散，本文取10。根据式(16)得到模型新的纵向网格间距后，需对初始模型沿深度方向非均匀采样，以保证模型总深度不变。采样过程如图 3所示，具体采样流程如下。

首先，从z=0开始给一个较小的试探步长(图 3中的第一个蓝色矩形)并不断增大该试探步长直至得到η₁深度处网格点，有η₁=Δz(η₁)；接着，从z=η₁开始重复上述过程，得到η₂深度处网格点，并有η₂－η₁=Δz(η₁)；最后，重复该过程直至达到模型最大深度，则得到一个纵向网格间距随地层速度自适应变化的新模型。新模型中网格点的速度与规则网格模型在同一深度下网格点的速度相同。此外，由于地下介质速度主要沿深度方向变化，而沿水平方向变化不大，因此，新模型的横向网格间距设为固定值。上述变网格方法是一种全局可变网格方法，能够实现网格大小随模型速度自适应变化，无须改变差分系数或在不同尺度的网格间设置过渡带。

根据本文采样方法，采样后的变网格模型可视作对初始模型进行坐标变换

$ \left\{\begin{array}{l} \xi=x \\ \eta=\varphi(z) \end{array}\right. $

(17)

式中：ξ和η为变换后坐标变量；φ(z)表示坐标系纵向映射关系，描述了图 3中所展示的自适应变网格采样方法。

应用变网格模型进行波场延拓时，原有的声波方程已不再适用，可利用映射关系推导与之相匹配的修正声波方程。空间坐标z关于η的一阶导数可写成

$ \frac{\partial z}{\partial \eta}=\frac{1}{\varphi^{\prime}(z)} $

(18)

二阶导数为

$ \frac{\partial^2 z}{\partial \eta^2}=\frac{\partial\left[\frac{1}{\varphi^{\prime}(z)}\right]}{\partial z} \frac{\partial z}{\partial \eta}=-\frac{\varphi^{\prime \prime}(z)}{\left[\varphi^{\prime}(z)\right]^3} $

(19)

根据链式法则，在式(5)中，声压u(x, z)关于空间坐标η的一阶导数可写成

$ \frac{\partial u}{\partial \eta}=\frac{\partial u}{\partial x} \frac{\partial x}{\partial \eta}+\frac{\partial u}{\partial z} \frac{\partial z}{\partial \eta}=\frac{\partial u}{\partial z} \frac{\partial z}{\partial \eta}=\frac{1}{\varphi^{\prime}(z)} \frac{\partial u}{\partial z} $

(20)

二阶导数可写成

$ \begin{aligned} & \frac{\partial^2 u}{\partial \eta^2}=\frac{\partial \frac{\partial u}{\partial \eta}}{\partial x} \frac{\partial x}{\partial \eta}+\frac{\partial \frac{\partial u}{\partial \eta}}{\partial z} \frac{\partial z}{\partial \eta} \\ & =\left(\frac{\partial^2 u}{\partial z^2} \frac{\partial z}{\partial \eta}+\frac{\partial \frac{\partial u}{\partial \eta}}{\partial z} \frac{\partial u}{\partial z}\right) \frac{\partial z}{\partial \eta} \\ & =\frac{\partial^2 u}{\partial z^2}\left(\frac{\partial z}{\partial \eta}\right)^2+\frac{\partial u}{\partial z} \frac{\partial^2 z}{\partial \eta^2} \\ & =\left[\frac{1}{\varphi^{\prime}(z)}\right]^2 \frac{\partial^2 u}{\partial z^2}-\frac{\varphi^{\prime \prime}(z)}{\left[\varphi^{\prime}(z)\right]^3} \frac{\partial u}{\partial z} \end{aligned} $

(21)

因此，式(5)可改写为

$ \frac{1}{v^2} \frac{\partial^2 u}{\partial t^2}=\frac{\partial^2 u}{\partial x^2}+\left[\frac{1}{\varphi^{\prime}(z)}\right]^2 \frac{\partial^2 u}{\partial z^2}-\frac{\varphi^{\prime \prime}(z)}{\left[\varphi^{\prime}(z)\right]^3} \frac{\partial u}{\partial z} $

(22)

与式(5)相比，式(22)需要额外计算声压关于空间坐标z的一阶导数，单个网格点的计算量有所增加，但由于φ″(z)远小于[φ′(z)]³，忽略该项几乎不影响计算精度，因而式(22)可进一步简化为

$ \frac{1}{v^2} \frac{\partial^2 u}{\partial t^2}=\frac{\partial^2 u}{\partial x^2}+\left[\frac{1}{\varphi^{\prime}(z)}\right]^2 \frac{\partial^2 u}{\partial z^2} $

(23)

式(23)即本文采用的变网格声波方程，不难发现，变网格模型仍需要在坐标系(x，z)中参与计算，这是因为坐标系(ξ，η)为不规则坐标系，难以在该坐标系中对波动方程进行编程求解。

式(6)可进一步改写为

$ \begin{aligned} u_{i, j}^{n+1}= & 2 u_{i, j}^n-u_{i, j}^{n-1}+\frac{v^2(\Delta t)^2}{(\Delta x)^2}\left(u_{*, j}^n\right)^{\prime \prime}+ \\ & \frac{v^2(\Delta t)^2}{(\Delta z)^2} \frac{1}{\mathit{\Psi }_j^2}\left(u_{i, *}^n\right)^{\prime \prime} \end{aligned} $

(24)

式中Ψ_j为φ′(z)的离散值，可采用有限差分法求解。

4 正演模拟 4.1 四层模型

为了验证本文方法的正确性及优势，设计了一个四层模型(图 4)。分别采用细网格、粗网格以及自适应变网格方法进行正演模拟测试。细网格模型如图 4a所示，细网格模型网格点数为401×401，横、纵向网格间距均为5 m；粗网格模型网格点数为401×251，横、纵向网格间距分别为5、8 m；变网格模型如图 4b所示，由细网格模型经自适应采样得到，网格点数为401×241，横向网格间距为5 m，纵向网格间距如图 5所示，沿深度方向自适应变化。由图 4b可以看出，经自适应采样后，模型低速浅层网格点数基本不变，而高速层被大幅压缩，有效避免了深层过采样现象。此外，从图 5可以看到，变网格模型的纵向网格间距能够在不同速度层之间平滑过渡，能够有效避免虚假反射现象。

在正演模拟中，声波方程的求解采用式(24)，震源采用主频为30 Hz的Ricker子波，位于模型表面中心，时间采样间隔为0.5 ms，采样时长为1.6 s。图 6是采用三种方案模拟的炮记录，如黑色箭头所示，细网格(图 6a)和变网格(图 6b)方法的合成记录基本没有出现数值频散，而粗网格炮记录(图 6c)数值频散较为明显(箭头所示)。在相同的计算资源条件下(计算机CPU型号为Intel(R) Xeon(R) Silver CPU 2.20 GHz，下同)，细网格、粗网格和变网格方法的计算时间分别为142、89和92 s。可见，采用本文自适应变网格方法既可以保证模拟精度，又可以提高计算效率、降低内存需求。

4.2 Marmousi模型

利用Marmousi模型测试本文方法对复杂模型适应性和计算效率。采用规则网格方法和自适应变网格方法进行正演模拟测试，规则网格模型如图 7a所示，网格点数为1360×700，横、纵向网格间距均为5 m；变网格模型如图 7b所示，由规则网格模型经自适应采样得到，横向网格点数和网格间距不变，纵向网格点数减少至436，减少了37.7%，纵向网格间距如图 8所示。可以看出，经自适应变网格采样后，Marmousi模型的高速层被压缩，网格点数大幅减少，网格间距随速度的增大而增大，但模型仍能保持原有的结构。

正演时的采样间隔为0.4 ms，震源采用主频为30 Hz的Ricker子波，布置在地表中心。图 9a和图 9b分别是采用规则网格和变网格方法得到的1.2 s时刻的模拟波场，图 9c展示了图 9b的线性插值结果，图 9d是图 9a与图 9c之差，相对于规则网格模拟波场，插值波场在L₂范数下的相对误差在1%以下。图 10a和图 10b分别展示了图 9a、图 9c和图 9d所示波场及波场差值在3.5 km和4.5 km处的单道曲线对比，不难看出，规则网格和变网格模拟波场的单道曲线存在着较小的振幅和相位误差，该误差主要源自式(16)中每个波长内的采样点数l的取值：增大l的值，可以降低模拟误差，但会削弱变网格方法加速计算的效果；减小l的值，可以进一步提高计算效率，但模拟波场的振幅和相位误差也会增大。在本例中，采用规则网格模型和变网格模型所需的计算时间分别为94、61 s，应用本文方法后效率提升了约35.1%。由此可见，本文自适应变网格方法仍适用于较复杂模型，能够降低内存占用、提高正演模拟效率。

5 逆时偏移成像

将本文基于自适应变网格策略的紧致差分数值模拟方法应用于Marmousi模型逆时偏移，以进一步验证本文方法的优越性。规则网格模型和变网格模型如图 7所示，成像参数如下：时间采样间隔为0.4 ms，总接收时长为3.0 s，震源为主频30 Hz的Ricker子波，采用全接收观测系统，地表布置68炮，炮间距为100 m，地表布置1360个接收点，接收点间距为5 m。采用一种震源归一化互相关成像条件对正、反传波场进行成像，成像条件为

$ I(x, z)=\frac{\int_0^{T_{\max }} S(x, z, t) R(x, z, t) \mathrm{d} t}{\sum\limits_0^{T_{\max }} S^2(x, z, t)} $

(25)

式中：I(x，z)为成像结果；S(x, z, t)为正向延拓震源波场；R(x, z, t)为反向延拓检波点波场；T_max为记录长度。

图 11a、图 11b分别展示了规则网格和自适应变网格方法的偏移剖面，图 11c为图 11b的线性插值结果。由图可见，两种方法均能成功地对复杂Marmousi模型进行成像。在计算效率方面，规则网格和变网格方法的计算时间分别为756、501 s，采用变网格后，计算效率提升了33.7%。抽取偏移结果在横向2、5 km处的单道曲线进行对比，如图 12所示，可以看到，两种方法的单道曲线仅有较小的振幅误差，基本没有深度误差。

综上所述，本文方法能够在保证Marmousi模型成像精度的前提下显著提高计算效率。

6 结论和讨论

本文构造了求解二维声波方程的时间二阶中心差分、空间四阶紧致差分离散格式，并引入自适应变网格策略进行计算加速，通过理论推导和模型测试，得到以下结论和认识：

(1) 空间四阶紧致差分格式的计算精度显著高于四阶中心差分格式，与六阶中心差分格式较为接近；

(2) 自适应变网格方法能够实现全局网格剖分优化，可在保证计算精度的同时提高计算效率、节约内存占用；

(3) 本文方法适用于复杂模型逆时偏移成像，与规则细网格成像结果误差较小。

此外，本文所采用的变网格方法仅考虑了深度方向，而横向上则仍使用规则网格，这是因为在常规情况下，地层速度主要沿深度方向变化，而且同时对横向和纵向进行变网格剖分将会带来一系列问题。首先，模型将不再是矩形，而是变得极为不规则，这使得网格剖分和离散计算较为困难；其次，地震波动方程可能增加若干额外项，变得更为复杂，导致单个网格点的计算量上升，反而起不到提高计算效率的效果；其三，模拟和成像精度可能会有所下降，影响最终结果。