一种基于隐式离散模型的频率域声波仿真方法及系统

1.本发明涉及声波技术领域，特别涉及一种基于隐式离散模型的频率域声波仿真方法及系统。


背景技术：

2.目前，随着超深检测复杂度的增加，对于高精度的成像方法的需求日益增大。由于全波形反演技术能利用超声记录中的各种波形信息来重建高精度的介质参数，成为国内外学者研究的热点问题。全波形反演主要分为频率域全波形反演，时间域全波形反演和混合域全波形反演；由于并行性好、多尺度明显以及能快速选频反演，频率域全波形反演被广泛应用研究。频率域仿真是频率域全波形反演的核心，研究高精度，高效的频率域仿真方法能提高反演的效率和反演的质量。
3.常规的频率域声波有限差分仿真方法，利用显式中心二阶有限差分方法来离散频率域声波方程。这种方法会产生严重的数值频散，并且仿真的结果精度很低。保持归一化相速度在1％的误差范围内，每个波长至少需要13个采样点，仿真的效率极低。为了提高仿真效率和仿真模拟的精度，以解决旋转坐标方法的局限性，chen(2012)提出一种平均导数的方法，该方法利用正交方向的加权平均有限差分算子来离散频率域声波方程，然后用matlab软件中的非线性优化器fmincon来计算优化系数。tang(2015)将平均导数用于高阶有限差分方法，利用模拟退火优化方法来优化系数，取得了一定的效果。借助于平均导数的思想，fan(2017)提出来了一种广义优化方法来近似拉普拉斯算子，并对不同的方法进行了对比。
4.综上所述，频率域声波有限差分仿真的改进方法，主要可分为旋转坐标法和平均导数法。旋转坐标法是利用旋转坐标算子和质量加速项来离散方程，得到优化差分系数。该类方法能在一定程度上降低频散，提高数值模拟的精度，但是它只适用于等间隔采样的情况，实际使用局限性很高。平均导数法是利用正交方向的加权平均有限差分算子来离散方程。该类方法能适用不等间隔采样的情况，并且对数值频散有一定的压制作用。具体地，基于旋转坐标法，通过组合不同的旋转坐标的差分算子来近似拉普拉斯算子，利用最小二乘优化方法求出优化系数，然后进行三维频率域声波方程仿真。基于平均导数法，利用正交方向的加权有限差分算子离散拉普拉斯算子，利用最优化方法求出优化系数，提高数值模拟的精度。
5.然而，上述现有的频率域声波仿真方法均借助于显式有限差分方法，即显示离散模型，精度有限，且数值频散严重；均需要优化大量的参数，仿真效率低；均需要结合质量加速技术才能提高数值模拟的精度，实际使用局限性高。


技术实现要素：

6.本发明的目的在于，克服现有基于显示离散模型的频率域声波仿真模拟精度低、数值频散严重和效率低的问题，从而提出一种一种基于隐式离散模型的频率域声波仿真方
法及系统，可以显著地压制数值频散，提高仿真模拟的精度，提高计算效率，节省仿真时间，为频率域全波形反演提供高精度、高效的仿真模拟方案，并且可以根据需要选择不同的方案，灵活方便，适用范围广。
7.为解决上述技术问题，本发明利用隐式离散模型来离散拉普拉斯算子，通过保留不同阶的结果，可以得到不同精度的离散方案，利用约束的信赖域方法获得差分系数，然后进行三维频率域声波方程仿真；相比于前人的方案，优化系数少，节省计算消耗，不需要质量加速技术，降低复杂度，可以根据需要选择不同的方案，灵活性更好，适用范围更广；本发明提供的一种基于隐式离散模型的频率域声波仿真方法，包括以下步骤：
8.步骤1)获取差分系数、声波波长的采样点数、采样间隔、声波传播角度和声波方位角；
9.步骤2)通过约束的信赖域算法，利用所述差分系数、所述声波波长的采样点数、所述采样间隔、所述声波传播角度和所述声波方位角获得最优差分系数；
10.步骤3)将所述最优差分系数输入预设的离散模型中，以获得最终离散模型；
11.步骤4)基于所述最终离散模型、频率域震源函数和吸收边界条件，构造阻抗矩阵；
12.步骤5)基于所述阻抗矩阵，获得频率域声波波场。
13.作为上述方法的一种改进，所述步骤2)具体包括以下迭代过程：
14.步骤2-1)基于所述差分系数、所述声波波长的采样点数、所述采样间隔、所述声波传播角度和所述声波方位角，确定归一化相速度误差的平方和；
15.步骤2-2)基于所述差分系数、所述声波波长的采样点数、所述声波传播角度和所述声波方位角，确定与所述归一化相速度误差的平方和对应的更新后的差分系数；
16.步骤2-3)判断当前迭代次数是否大于预设迭代次数；如果不是，则将所述更新后的差分系数作为下次迭代中的差分系数，令迭代次数加一，返回所述步骤2-1)；如果是，执行步骤2-4)；
17.步骤2-4)从多次迭代获取的所述归一化相速度误差的平方和中选取最小值，将所述最小值对应的所述更新后的差分系数作为最优差分系数。
18.作为上述方法的一种改进，确定归一化相速度误差的平方和包括：
19.根据所述差分系数、所述声波波长的采样点数、所述采样间隔、所述声波传播角度和所述声波方位角确定归一化相速度；
20.根据所述归一化相速度确定归一化相速度误差的平方和。
21.作为上述方法的一种改进，所述步骤3)具体包括：
22.基于所述声波波长的采样点数，从若干初始离散模型中选择最优离散模型；
23.将所述最优差分系数输入预设的最优离散模型中，以获得最终离散模型。
24.作为上述方法的一种改进，所述若干初始离散模型包括：隐式7点离散模型、隐式19点离散模型和隐式27点离散模型。
25.作为上述方法的一种改进，所述步骤5)具体包括：
26.采用多波前lu分解的方式分解所述阻抗矩阵，得到频率域声波波场。
27.为实现本发明再一目的，本发明一种基于隐式离散模型的频率域声波仿真系统，用于执行上述基于隐式离散模型的频率域声波仿真方法；所述系统包括：获取单元、迭代单元、最终离散模型单元、阻抗矩阵构造单元和声波波场单元；其中，
28.所述获取单元，用于获取差分系数、声波波长的采样点数、采样间隔、声波传播角度和声波方位角；
29.所述迭代单元，用于执行约束的信赖域算法，利用所述差分系数、所述声波波长的采样点数、所述采样间隔、所述声波传播角度和所述声波方位角获得最优差分系数；
30.所述最终离散模型单元，用于将所述最优差分系数输入预设的离散模型中，得到最终离散模型；
31.所述阻抗矩阵构造单元，用于根据所述最终离散模型、频率域震源函数和吸收边界条件构造阻抗矩阵；
32.所述声波波场单元，用于根据所述阻抗矩阵得到频率域声波波场。
33.作为上述系统的一种改进，所述迭代单元用于执行如下迭代处理：
34.基于所述差分系数、所述声波波长的采样点数、所述采样间隔、所述声波传播角度和所述声波方位角，用于确定归一化相速度误差的平方和；
35.基于所述差分系数、所述声波波长的采样点数、所述声波传播角度和所述声波方位角，用于确定与所述归一化相速度误差的平方和对应的更新后的差分系数；
36.判断当前迭代次数是否大于预设迭代次数；
37.如果不是，则利用所述更新后的差分系数作为下次迭代中的差分系数，并重新进行迭代处理，同时令迭代次数加一；
38.如果是，则从多次迭代获取的所述归一化相速度误差的平方和中选取最小值，将所述最小值对应的所述更新后的差分系数作为最优差分系数。
39.作为上述系统的一种改进，所述迭代单元用于：
40.根据所述差分系数、所述声波波长的采样点数、所述采样间隔、所述声波传播角度和所述声波方位角确定归一化相速度；
41.根据所述归一化相速度确定归一化相速度误差的平方和。
42.作为上述系统的一种改进，所述系统还包括：最优离散模型单元，其中，
43.所述最优离散模型单元，用于根据所述声波波长的采样点数从若干初始离散模型中选择最优离散模型；
44.所述最终离散模型单元具体用于：
45.将所述最优差分系数输入预设的最优离散模型中，得到最终离散模型。
46.作为上述系统的一种改进，所述声波波场单元具体用于：
47.采用多波前lu分解的方式分解所述阻抗矩阵，得到频率域声波波场。
48.本发明还提供一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现上述基于隐式离散模型的频率域声波仿真方法的步骤。
49.本发明还提供一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现上述基于隐式离散模型的频率域声波仿真方法的步骤。
50.本发明采用了隐式离散模型，隐式离散模型的结果与理论离散模型吻合度高，精度高，可以有效压制三维频率域声波正演模拟的数值频散，提高三维频率域声波仿真的精度，并且随着保留阶数的增高，隐式离散模型压制数值频散的效果更好。并且采用隐式离散模型还可以减少优化参数，节省资源，因此，扩大了适用范围。
附图说明
51.图1是本发明实施例中基于隐式离散模型的频率域声波仿真方法的流程图；
52.图2是本发明实施例中隐式7点离散模型的格式图；
53.图3是本发明实施例中隐式19点离散模型的格式图；
54.图4是本发明实施例中隐式27点离散模型的格式图；
55.图5是η＝30
°
时现有7点离散模型的频散曲线图；
56.图6是本发明实施例中η＝30
°
时隐式7点离散模型的频散曲线图；
57.图7是本发明实施例中η＝30
°
时隐式19离散模型的频散曲线图；
58.图8是本发明实施例中η＝30
°
时隐式27点离散模型的频散曲线图；
59.图9是θ＝30
°
时现有7点离散模型的频散曲线图；
60.图10是本发明实施例中θ＝30
°
时隐式7点离散模型的频散曲线图；
61.图11是本发明实施例中θ＝30
°
时隐式19离散模型的频散曲线图；
62.图12是本发明实施例中θ＝30
°
时隐式27点离散模型的频散曲线图；
63.图13是现有7点离散模型在0.1s时刻的时间域波场快照图；
64.图14是本发明实施例中隐式7点离散模型在0.1s时刻的时间域波场快照图；
65.图15是本发明实施例中隐式19点离散模型在0.1s时刻的时间域波场快照图；
66.图16是本发明实施例中隐式27点离散模型在0.1s时刻的时间域波场快照图；
67.图17是理论离散模型的时间域地震记录图；
68.图18是现有7点离散模型的时间域地震记录图；
69.图19是本发明实施例中隐式7点离散模型的时间域地震记录图；
70.图20是本发明实施例中隐式19点离散模型的时间域地震记录图；
71.图21是本发明实施例中隐式27点离散模型的时间域地震记录图；
72.图22是本发明实施例中基于隐式离散模型的频率域声波仿真系统的结构框图。
具体实施方式
73.下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
74.本领域技术人员知道，本发明的实施方式可以实现为一种系统、装置、设备、方法或计算机程序产品。因此，本公开可以具体实现为以下形式，即：完全的硬件、完全的软件(包括固件、驻留软件、微代码等)，或者硬件和软件结合的形式。
75.鉴于现有技术耗费了大量资源，实际使用的局限性很高且精度有限，本发明实施例提供一种基于隐式离散模型的频率域声波仿真方法，以节省资源，扩大适用范围，增加精度。以下结合附图对本发明进行详细说明。
76.图1是本发明实施例中基于隐式离散模型的频率域声波仿真方法的流程图。如图1所示，基于隐式离散模型的频率域声波仿真方法包括：
77.s101：获取差分系数、声波波长的采样点数、采样间隔、声波传播角度和声波方位角。
78.执行如下迭代处理：
79.步骤1)获取差分系数、声波波长的采样点数、采样间隔、声波传播角度和声波方位角；
80.步骤2)通过约束的信赖域算法，利用所述差分系数、所述声波波长的采样点数、所述采样间隔、所述声波传播角度和所述声波方位角获得最优差分系数；
81.所述步骤2)具体包括以下迭代过程：
82.步骤2-1)基于所述差分系数、所述声波波长的采样点数、所述采样间隔、所述声波传播角度和所述声波方位角，确定归一化相速度误差的平方和；
83.步骤2-2)基于所述差分系数、所述声波波长的采样点数、所述声波传播角度和所述声波方位角，确定与所述归一化相速度误差的平方和对应的更新后的差分系数；其中，确定归一化相速度误差的平方的步骤包括：
84.根据所述差分系数、所述声波波长的采样点数、所述采样间隔、所述声波传播角度和所述声波方位角确定归一化相速度；
85.根据所述归一化相速度确定归一化相速度误差的平方和；
86.步骤2-3)判断当前迭代次数是否大于预设迭代次数；如果不是，则将所述更新后的差分系数作为下次迭代中的差分系数，令迭代次数加一，返回所述步骤2-1)；如果是，执行步骤2-4)；
87.步骤2-4)从多次迭代获取的所述归一化相速度误差的平方和中选取最小值，将所述最小值对应的所述更新后的差分系数作为最优差分系数。
88.步骤3)将所述最优差分系数输入预设的离散模型中，以获得最终离散模型；
89.步骤4)基于所述最终离散模型、频率域震源函数和吸收边界条件，构造阻抗矩阵；
90.步骤5)基于所述阻抗矩阵，获得频率域声波波场。
91.图1所示的基于隐式离散模型的频率域声波仿真方法的执行主体可以为计算机。由图1所示的流程可知，本发明实施例的基于隐式离散模型的频率域声波仿真方法先根据差分系数确定归一化相速度误差的平方和更新后的差分系数；当迭代次数大于预设迭代次数时，从多个归一化相速度误差的平方和中选取最小值，将最小值对应的差分系数作为最优差分系数，否则将更新后的差分系数作为下次迭代中的差分系数，令迭代次数加一，重新执行迭代处理；然后将最优差分系数输入预设的离散模型中，得到最终离散模型；接着根据最终离散模型、频率域震源函数和吸收边界条件构造阻抗矩阵，最后根据阻抗矩阵得到频率域声波波场，可以减少优化参数，节省资源，扩大适用范围，增加精度，进一步提高了声波正演的效率和质量。
92.一实施例中，步骤2-2)包括：
93.根据差分系数、声波波长的采样点数、采样间隔、声波传播角度和声波方位角确定归一化相速度。
94.其中，通过如下公式确定归一化相速度：
[0095][0096]
a＝(2e
x-2)+βy(2e
x-2)(2e
y-2)+βz(2e
x-2)(2e
z-2)+βyβz(2e
x-2)(2e
y-2)(2e
z-2)；
[0097]
b＝(2e
y-2)+β
x
(2e
y-2)(2e
x-2)+βz(2e
y-2)(2e
z-2)+βzβ
x
(2e
x-2)(2e
y-2)(2e
z-2)；
[0098]
c＝(2e
z-2)+β
x
(2e
z-2)(2e
x-2)+βy(2e
z-2)(2e
y-2)+β
x
βy(2e
x-2)(2e
y-2)(2e
z-2)；
[0099]
d＝(1+β
x
(2e
x-2)+βy(2e
y-2)+βz(2e
z-2)+β
x
βy(2e
x-2)(2e
y-2)+
[0100]
βyβz(2e
y-2)(2e
z-2)+βzβ
x
(2e
z-2)(2e
x-2)+β
x
βyβz(2e
x-2)(2e
y-2)(2e
z-2))；
[0101][0102][0103][0104][0105]
其中，v
ph
为相速度，v为声波在介质中的传播速度，为归一化相速度，g为声波波长的采样点数，β
x
为x轴方向的差分系数，βy为y轴方向的差分系数，βz为z轴方向的差分系数，δx为x轴方向的采样间隔，δy为y轴方向的采样间隔，δz为z轴方向的采样间隔，θ为声波传播角度，η为声波方位角。θ和η的取值范围均是[0,π/2]。
[0106]
根据归一化相速度确定归一化相速度误差的平方和。
[0107]
其中，通过如下公式确定归一化相速度误差的平方和：
[0108][0109]
其中，e为归一化相速度误差的平方和。
[0110]
一实施例中，在执行步骤3)之前，还包括：根据声波波长的采样点数从多个初始离散模型中选择最优离散模型。
[0111]
步骤3)包括：将最优差分系数输入预设的最优离散模型中，得到最终离散模型。
[0112]
其中，初始离散模型包括隐式7点离散模型、隐式19点离散模型和隐式27点离散模型。
[0113]
隐式7点离散模型如下：
[0114][0115]
其中，δx为x轴方向的采样间隔，δy为y轴方向的采样间隔，δz为z轴方向的采样间隔，β
x
为x轴方向的差分系数，βy为y轴方向的差分系数，βz为z轴方向的差分系数，ω为圆频率，v为波在介质中的传播速度，p
i+1,j,k
为离散点(i+1,j,k)的波场，p
i-1,j,k
为离散点(i-1,j,k)的波场，p
i,j+1,k
为离散点(i,j+1,k)的波场，p
i,j-1,k
为离散点(i,j-1,k)的波场，p
i,j,k+1
为离散点(i,j,k+1)的波场，p
i,j,k-1
为离散点(i,j,k-1)的波场，p
i,j,k
为离散点(i,j,k)的波
场。
[0116]
图2是本发明实施例中隐式7点离散模型的格式图。如图2所示，隐式7点离散模型只保留了二阶的有限差分算子，只用到立方体的中心点和六个面内的点。
[0117]
隐式19点离散模型如下：
[0118][0119]
其中，p
i+1,j+1,k
为离散点(i+1,j+1,k)的波场，p
i+1,j-1,k
为离散点(i+1,j-1,k)的波场，p
i-1,j+1,k
为离散点(i-1,j+1,k)的波场，p
i-1,j-1,k
为离散点(i-1,j-1,k)的波场，p
i,j+1,k+1
为离散点(i,j+1,k+1)的波场，p
i,j+1,k-1
为离散点(i,j+1,k-1)的波场，p
i,j-1,k+1
为离散点(i,j-1,k+1)的波场，p
i,j-1,k-1
为离散点(i,j-1,k-1)的波场，p
i+1,j,k+1
为离散点(i+1,j,k+1)的波场，p
i-1,j,k+1
为离散点(i-1,j,k+1)的波场，p
i+1,j,k-1
为离散点(i+1,j,k-1)的波场，p
i-1,j,k-1
为离散点(i-1,j,k-1)的波场。
[0120]
图3是本发明实施例中隐式19点离散模型的格式图。如图3所示，隐式19点离散模型只保留了二阶和四阶的有限差分算子，它不但用到了立方体的中心点和六个面内的点，还用到12条棱上的点，这是计算四阶有限差分算子所需要的额外的点。
[0121]
隐式27点离散模型如下：
[0122][0123]
其中，p
i+1,j+1,k+1
为离散点(i+1,j+1,k+1)的波场，p
i+1,j-1,k+1
为离散点(i+1,j-1,k+1)的波场，p
i-1,j+1,k+1
为离散点(i-1,j+1,k+1)的波场，p
i-1,j-1,k+1
为离散点(i-1,j-1,k+1)的波场，p
i+1,j+1,k-1
为离散点(i+1,j+1,k-1)的波场，p
i+1,j-1,k-1
为离散点(i+1,j-1,k-1)的波场，p
i-1,j+1,k-1
为离散点(i-1,j+1,k-1)的波场，p
i-1,j-1,k-1
为离散点(i-1,j-1,k-1)的波场。
[0124]
图4是本发明实施例中隐式27点离散模型的格式图。如图4所示，隐式27点离散模型保留了二阶、四阶和六阶的有限差分算子，它不但用到了立方体的中心点和六个面内的点，还用到12条棱上的点，这是计算六阶有限差分算子所需要的额外的点。
[0125]
本发明实施例的具体流程如下：
[0126]
1、根据差分系数、声波波长的采样点数、采样间隔、声波传播角度和声波方位角确定归一化相速度。
[0127]
2、根据归一化相速度确定归一化相速度误差的平方和。
[0128]
3、根据差分系数、声波波长的采样点数、声波传播角度和声波方位角确定更新后
的差分系数。
[0129]
4、判断迭代次数是否大于预设迭代次数；当大于预设迭代次数时，从多个归一化相速度误差的平方和中选取最小值，将最小值对应的差分系数作为最优差分系数；当小于或等于预设迭代次数时，将更新后的差分系数作为下次迭代中的差分系数，令迭代次数加一，返回步骤1。
[0130]
5、根据声波波长的采样点数从隐式7点离散模型、隐式19点离散模型和隐式27点离散模型中选择最优离散模型。
[0131]
表1
[0132][0133]
表2
[0134][0135]
表3
[0136][0137]
表1是第一实施例中隐式7点离散模型对应的最优差分系数表。表2第二实施例中是隐式19点离散模型对应的最优差分系数表。表3是第三实施例中隐式27点离散模型对应的最优差分系数表。
[0138]
图5是η＝30
°
时现有7点离散模型的频散曲线图。图6是本发明实施例中η＝30
°
时隐式7点离散模型的频散曲线图。图7是本发明实施例中η＝30
°
时隐式19离散模型的频散曲线图。图8是本发明实施例中η＝30
°
时隐式27点离散模型的频散曲线图。图5-图8的横坐标均为采样点数的倒数纵坐标均为归一化相速度图中的四条频散曲线分别为θ等于0
°
、30
°
、60
°
和90
°
时的频散曲线。
[0139]
图9是θ＝30
°
时现有7点离散模型的频散曲线图。图10是本发明实施例中θ＝30
°
时隐式7点离散模型的频散曲线图。图11是本发明实施例中θ＝30
°
时隐式19离散模型的频散曲线图。图12是本发明实施例中θ＝30
°
时隐式27点离散模型的频散曲线图。图9-图12的横坐标均为采样点数的倒数纵坐标均为归一化相速度图中的四条频散曲线分别为η等于0
°
、30
°
、60
°
和90
°
时的频散曲线。
[0140]
从图5-图12中可以看出，相较于现有7点离散模型，隐式离散模型可以显著地降低数值频散，提高正演模拟的精度。随着保留阶数的增高，隐式离散模型压制数值频散的效果更好。
[0141]
6、将最优差分系数输入预设的最优离散模型中，得到最终离散模型。
[0142]
7、根据最终离散模型、频率域震源函数和吸收边界条件构造阻抗矩阵。
[0143]
8、采用多波前lu分解的方式分解阻抗矩阵，得到频率域声波波场。
[0144]
在其中一种实施例中，设波场模型的v＝3500m/s，x轴方向有50个网格，y轴方向有50个网格，z轴方向有40个网格，三个方向的采样间隔均为15m，波场模型周围有10层pml吸收边界层，频率域震源函数是主频为30hz的ricker子波。图13是现有7点离散模型在0.1s时刻的时间域波场快照图。图14是本发明实施例中隐式7点离散模型在0.1s时刻的时间域波场快照图。图15是本发明实施例中隐式19点离散模型在0.1s时刻的时间域波场快照图。图16是本发明实施例中隐式27点离散模型在0.1s时刻的时间域波场快照图。从图13-图16中
可以看出，现有7点离散模型的频散严重，隐式离散模型可以显著地降低数值频散。随着保留阶数的增高，隐式离散模型压制数值频散的效果更好。图17是理论离散模型的时间域地震记录图。图18是现有7点离散模型的时间域地震记录图。图19是本发明实施例中隐式7点离散模型的时间域地震记录图。图20是本发明实施例中隐式19点离散模型的时间域地震记录图。图21是本发明实施例中隐式27点离散模型的时间域地震记录图。从图17-图21中可以看出，相较于现有7点离散模型，隐式离散模型的结果与理论离散模型吻合得更好，说明隐式离散模型的精度很高，尤其是隐式27点离散模型与理论离散模型吻合度非常高。这些结果证明了本技术可以有效压制三维频率域声波正演模拟的数值频散，提高三维频率域声波正演模拟的精度。
[0145]
综上，本发明实施例的基于隐式离散模型的频率域声波仿真方法先根据差分系数确定归一化相速度误差的平方和更新后的差分系数；当迭代次数大于预设迭代次数时，从多个归一化相速度误差的平方和中选取最小值，将最小值对应的差分系数作为最优差分系数，否则将更新后的差分系数作为下次迭代中的差分系数，令迭代次数加一，重新执行迭代处理；然后将最优差分系数输入预设的离散模型中，得到最终离散模型；接着根据最终离散模型、频率域震源函数和吸收边界条件构造阻抗矩阵，最后根据阻抗矩阵得到频率域声波波场，可以减少优化参数，节省资源，扩大适用范围，增加精度，进一步提高了声波正演的效率和质量。
[0146]
基于同一发明构思，本发明实施例还提供了一种基于隐式离散模型的频率域声波仿真系统，由于该系统解决问题的原理与基于隐式离散模型的频率域声波仿真方法相似，因此该系统的实施可以参见方法的实施，重复之处不再赘述。
[0147]
图22是本发明实施例中基于隐式离散模型的频率域声波仿真系统的结构框图。如图22所示，基于隐式离散模型的频率域声波仿真系统包括：
[0148]
获取单元，用于获取差分系数、声波波长的采样点数、采样间隔、声波传播角度和声波方位角；
[0149]
迭代单元，用于执行如下迭代处理：
[0150]
根据差分系数、声波波长的采样点数、采样间隔、声波传播角度和声波方位角确定归一化相速度误差的平方和；
[0151]
根据差分系数、声波波长的采样点数、声波传播角度和声波方位角确定更新后的差分系数；
[0152]
基于所述差分系数、所述声波波长的采样点数、所述采样间隔、所述声波传播角度和所述声波方位角，用于确定归一化相速度误差的平方和；
[0153]
基于所述差分系数、所述声波波长的采样点数、所述声波传播角度和所述声波方位角，用于确定与所述归一化相速度误差的平方和对应的更新后的差分系数；
[0154]
判断当前迭代次数是否大于预设迭代次数；
[0155]
如果不是，则利用所述更新后的差分系数作为下次迭代中的差分系数，并重新进行迭代处理，同时令迭代次数加一；
[0156]
如果是，则从多次迭代获取的所述归一化相速度误差的平方和中选取最小值，将所述最小值对应的所述更新后的差分系数作为最优差分系数。
[0157]
最终离散模型单元，用于将最优差分系数输入预设的离散模型中，得到最终离散
模型；
[0158]
阻抗矩阵构造单元，用于根据最终离散模型、频率域震源函数和吸收边界条件构造阻抗矩阵；
[0159]
声波波场单元，用于根据阻抗矩阵得到频率域声波波场。
[0160]
在其中一种实施例中，还包括：
[0161]
最优离散模型单元，用于根据声波波长的采样点数从多个初始离散模型中选择最优离散模型；
[0162]
最终离散模型单元具体用于：
[0163]
将最优差分系数输入预设的最优离散模型中，得到最终离散模型。
[0164]
在其中一种实施例中，迭代单元具体用于：
[0165]
根据差分系数、声波波长的采样点数、采样间隔、声波传播角度和声波方位角确定归一化相速度；
[0166]
根据归一化相速度确定归一化相速度误差的平方和。
[0167]
在其中一种实施例中，声波波场单元具体用于：
[0168]
采用多波前lu分解的方式分解阻抗矩阵，得到频率域声波波场。
[0169]
综上，本发明实施例的基于隐式离散模型的频率域声波仿真系统先根据差分系数确定归一化相速度误差的平方和更新后的差分系数；当迭代次数大于预设迭代次数时，从多个归一化相速度误差的平方和中选取最小值，将最小值对应的差分系数作为最优差分系数，否则将更新后的差分系数作为下次迭代中的差分系数，令迭代次数加一，重新执行迭代处理；然后将最优差分系数输入预设的离散模型中，得到最终离散模型；接着根据最终离散模型、频率域震源函数和吸收边界条件构造阻抗矩阵，最后根据阻抗矩阵得到频率域声波波场，可以减少优化参数，节省资源，扩大适用范围，增加精度，进一步提高了声波正演的效率和质量。
[0170]
本发明实施例还提供了一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，处理器执行计算机程序时可以实现基于隐式离散模型的频率域声波仿真方法的全部或部分内容，例如，处理器执行计算机程序时可以实现如下内容：
[0171]
获取差分系数、声波波长的采样点数、采样间隔、声波传播角度和声波方位角；
[0172]
执行如下迭代处理：
[0173]
根据差分系数、声波波长的采样点数、采样间隔、声波传播角度和声波方位角确定归一化相速度误差的平方和；
[0174]
根据差分系数、声波波长的采样点数、声波传播角度和声波方位角确定更新后的差分系数；
[0175]
基于所述差分系数、所述声波波长的采样点数、所述采样间隔、所述声波传播角度和所述声波方位角，用于确定归一化相速度误差的平方和；
[0176]
基于所述差分系数、所述声波波长的采样点数、所述声波传播角度和所述声波方位角，用于确定与所述归一化相速度误差的平方和对应的更新后的差分系数；
[0177]
判断当前迭代次数是否大于预设迭代次数；
[0178]
如果不是，则利用所述更新后的差分系数作为下次迭代中的差分系数，并重新进
行迭代处理，同时令迭代次数加一；
[0179]
如果是，则从多次迭代获取的所述归一化相速度误差的平方和中选取最小值，将所述最小值对应的所述更新后的差分系数作为最优差分系数；将最优差分系数输入预设的离散模型中，得到最终离散模型；
[0180]
根据最终离散模型、频率域震源函数和吸收边界条件构造阻抗矩阵；
[0181]
根据阻抗矩阵得到频率域声波波场。
[0182]
综上，本发明实施例的计算机设备先根据差分系数确定归一化相速度误差的平方和更新后的差分系数；当迭代次数大于预设迭代次数时，从多个归一化相速度误差的平方和中选取最小值，将最小值对应的差分系数作为最优差分系数，否则将更新后的差分系数作为下次迭代中的差分系数，令迭代次数加一，重新执行迭代处理；然后将最优差分系数输入预设的离散模型中，得到最终离散模型；接着根据最终离散模型、频率域震源函数和吸收边界条件构造阻抗矩阵，最后根据阻抗矩阵得到频率域声波波场，可以减少优化参数，节省资源，扩大适用范围，增加精度，进一步提高了声波正演的效率和质量。
[0183]
本发明实施例还提供了一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执行时可以实现基于隐式离散模型的频率域声波仿真方法的全部或部分内容，例如，处理器执行计算机程序时可以实现如下内容：
[0184]
获取差分系数、声波波长的采样点数、采样间隔、声波传播角度和声波方位角；
[0185]
执行如下迭代处理：
[0186]
根据差分系数、声波波长的采样点数、采样间隔、声波传播角度和声波方位角确定归一化相速度误差的平方和；
[0187]
根据差分系数、声波波长的采样点数、声波传播角度和声波方位角确定更新后的差分系数；
[0188]
基于所述差分系数、所述声波波长的采样点数、所述采样间隔、所述声波传播角度和所述声波方位角，用于确定归一化相速度误差的平方和；
[0189]
基于所述差分系数、所述声波波长的采样点数、所述声波传播角度和所述声波方位角，用于确定与所述归一化相速度误差的平方和对应的更新后的差分系数；
[0190]
判断当前迭代次数是否大于预设迭代次数；
[0191]
如果不是，则利用所述更新后的差分系数作为下次迭代中的差分系数，并重新进行迭代处理，同时令迭代次数加一；
[0192]
如果是，则从多次迭代获取的所述归一化相速度误差的平方和中选取最小值，将所述最小值对应的所述更新后的差分系数作为最优差分系数。
[0193]
将最优差分系数输入预设的离散模型中，得到最终离散模型；
[0194]
根据最终离散模型、频率域震源函数和吸收边界条件构造阻抗矩阵；
[0195]
根据阻抗矩阵得到频率域声波波场。
[0196]
综上，本发明实施例的计算机可读存储介质先根据差分系数确定归一化相速度误差的平方和更新后的差分系数；当迭代次数大于预设迭代次数时，从多个归一化相速度误差的平方和中选取最小值，将最小值对应的差分系数作为最优差分系数，否则将更新后的差分系数作为下次迭代中的差分系数，令迭代次数加一，重新执行迭代处理；然后将最优差分系数输入预设的离散模型中，得到最终离散模型；接着根据最终离散模型、频率域震源函
数和吸收边界条件构造阻抗矩阵，最后根据阻抗矩阵得到频率域声波波场，可以减少优化参数，节省资源，扩大适用范围，增加精度，进一步提高了声波正演的效率和质量。
[0197]
以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
[0198]
本领域技术人员还可以了解到本发明实施例列出的各种说明性逻辑块(illustrative logical block)，单元，和步骤可以通过电子硬件、电脑软件，或两者的结合进行实现。为清楚展示硬件和软件的可替换性(interchangeability)，上述的各种说明性部件(illustrative components)，单元和步骤已经通用地描述了它们的功能。这样的功能是通过硬件还是软件来实现取决于特定的应用和整个系统的设计要求。本领域技术人员可以对于每种特定的应用，可以使用各种方法实现所述的功能，但这种实现不应被理解为超出本发明实施例保护的范围。
[0199]
本发明实施例中所描述的各种说明性的逻辑块，或单元，或装置都可以通过通用处理器，数字信号处理器，专用集成电路(asic)，现场可编程门阵列或其它可编程逻辑装置，离散门或晶体管逻辑，离散硬件部件，或上述任何组合的设计来实现或操作所描述的功能。通用处理器可以为微处理器，可选地，该通用处理器也可以为任何传统的处理器、控制器、微控制器或状态机。处理器也可以通过计算装置的组合来实现，例如数字信号处理器和微处理器，多个微处理器，一个或多个微处理器联合一个数字信号处理器核，或任何其它类似的配置来实现。
[0200]
本发明实施例中所描述的方法或算法的步骤可以直接嵌入硬件、处理器执行的软件模块、或者这两者的结合。软件模块可以存储于ram存储器、闪存、rom存储器、eprom存储器、eeprom存储器、寄存器、硬盘、可移动磁盘、cd-rom或本领域中其它任意形式的存储媒介中。示例性地，存储媒介可以与处理器连接，以使得处理器可以从存储媒介中读取信息，并可以向存储媒介存写信息。可选地，存储媒介还可以集成到处理器中。处理器和存储媒介可以设置于asic中，asic可以设置于用户终端中。可选地，处理器和存储媒介也可以设置于用户终端中的不同的部件中。
[0201]
在一个或多个示例性的设计中，本发明实施例所描述的上述功能可以在硬件、软件、固件或这三者的任意组合来实现。如果在软件中实现，这些功能可以存储与电脑可读的媒介上，或以一个或多个指令或代码形式传输于电脑可读的媒介上。电脑可读媒介包括电脑存储媒介和便于使得让电脑程序从一个地方转移到其它地方的通信媒介。存储媒介可以是任何通用或特殊电脑可以接入访问的可用媒体。例如，这样的电脑可读媒体可以包括但不限于ram、rom、eeprom、cd-rom或其它光盘存储、磁盘存储或其它磁性存储装置，或其它任何可以用于承载或存储以指令或数据结构和其它可被通用或特殊电脑、或通用或特殊处理器读取形式的程序代码的媒介。此外，任何连接都可以被适当地定义为电脑可读媒介，例如，如果软件是从一个网站站点、服务器或其它远程资源通过一个同轴电缆、光纤电缆、双绞线、数字用户线(dsl)或以例如红外、无线和微波等无线方式传输的也被包含在所定义的电脑可读媒介中。所述的碟片(disk)和磁盘(disc)包括压缩磁盘、镭射盘、光盘、dvd、软盘和蓝光光盘，磁盘通常以磁性复制数据，而碟片通常以激光进行光学复制数据。上述的组合
也可以包含在电脑可读媒介中。
[0202]
最后所应说明的是，以上实施例仅用以说明本发明的技术方案而非限制。尽管参照实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，对本发明的技术方案进行修改或者等同替换，都不脱离本发明技术方案的精神和范围，其均应涵盖在本发明的权利要求范围当中。