基于POD有限元法的粘滞声波方程震源波场重构

宋淳 王璐 冯民富

叠前逆时偏移是地震勘探中一种流行的地下结构成像方法，其成像条件需要同时刻的震源波场值与检波器波场值. 这在实际计算中就需要把正演模拟的所有时刻的震源波场数据全部存储下来，存储需求大. 虽然震源波场重构技术可以降低对于波场数据的存储需求，但会引入额外的计算复杂度. 为解决这个问题，本文提出了POD有限元法，并将其应用于粘滞震源波场重构.这里的本征正交分解（Proper Orthogonal Decomposition， POD）方法是一种降维方法，能够在降低数据量的同时提供足够的计算精度.   数值算例显示，该方法比传统的有限元法更节省存储空间，能够加快重构速度.

震源波场重构； 粘滞声波方程； 有限元法； 本征正交分解

O241.82A2023.031005

收稿日期： 2022-05-18

基金项目： 国家自然科学基金（11971337）

作者简介： 宋淳 （2000-）， 男， 硕士研究生， 主要研究方向为计算数学. E-mail： 2234878330@qq.com

通讯作者： 王璐. E-mail： 530441397@qq.com

POD finite element method for source wave field reconstruction of viscous acoustic wave equations

SONG Chun1， WANG Lu2， FENG Min-Fu2

（1. School of Biomass Science and Engineering， Sichuan University， Chengdu 610064， China;

2. School of Mathematics， Sichuan University， Chengdu 610064， China）

Prestack reverse time migration is popular for imaging underground structures in seismic exploration. Its imaging condition requires to gather data from the source wavefield and receiver wavefield simultaneously， which means that we have to store all source wavefield data at all times of the forward simulation and huge storage demand. The source wavefield reconstruction technologies can be used to solve this problem at the cost of introducing additional computational complexity. In this paper we introduce the POD finite element method and apply it to the reconstruction of viscous source wavefield， here the Proper Orthogonal Decomposition （POD） method is used to reduce the storage demand by decreasing the dimensionality of data while keeping high enough accuracy. Numerical examples show that， compared with the traditional finite element method， the introduced method can save more computer memory and speed up the reconstruction.

Reconstruction of source wavefield; Viscous acoustic wave equation; Finite element method; Proper orthogonal decomposition method

（2010 MSC 65M60）

1 引 言

叠前逆时偏移方法［1，2］是一种流行的复杂地质体成像手段. 该方法可分为三部分，即震源波场的正向延拓，检波器波场的逆时延拓以及成像条件. 其中，成像条件一般采用互相关成像，因而需要震源波场与检波器波场同一时刻的波场值. 但是，由于震源波场外推和检波器波场外推在时间上并不同步，往往需要将所有时间点的震源波场都存储下来才能成像，进而导致计算机存储量剧增. 为解决这个问题，Dussaud等［3］首先提出了震源波场重构法. 该方法先正向外推一次震源波场，然后再进行一次震源波场的逆时外推，以保证震源波场和检波器波场外推在时间上同步，减少存储开销.

经过不断地改进，目前重构震源波场方法已有多种［4］. 常见的重构方法基于有限差分［5］. 有限差分法适合处理规则区域问题，且编程简单，因而得到广泛应用，但它不能处理复杂地质结构及复杂边界条件. 与有限差分法不同，有限元法［6，7］具有网格剖分灵活、对复杂结构区域容易求解等优势，且基于变分原理时该方法还可以处理复杂边界条件，在重构震源波场时有潜在的优势. 其中，有限元的边界值法是重构震源波场的常用方法［8］. 该方法在震源波场正演过程中只需存储最后两个时间层的波场值作为震源波场重构的初值. 然而，经典的有限元法计算量大. 随着网格剖分加密、计算时间增加，该方法需要求解大规模线性方程组，从而引入额外的计算复杂度，这个缺点就愈加明显.

POD方法是一种能够提供足够精度而计算自由度较少的降维方法. 通过构造POD基和低维POD空间，基于POD的数值方法的计算自由度远远小于有限元方法，因而能够简化计算、节省计算时间. 例如，Luo等［9］将POD方法推广应用于偏微分方程的数值计算格式，包括基于POD降维的有限差分法、有限元法、有限体积法等. 大量研究表明，将POD方法应用于经典有限元法会在保持有限元法优势的同时大大节省计算时间.

大多数叠前逆时偏移方法基于声波方程，该方程模拟了理想情况下的地震波传播规律. 然而，真实的地下介质具有粘滞性，常规的声波偏移方法并不能正确地处理地下结构成像，因而研究粘滞声波方程的数值解法是有意义的. 值得注意的是，能够表达粘滞性的波动方程有多种［10］，其中Deng和McMechan［11］在声波方程中添加波场的一阶时间导数来调整振幅，使方程在波场正演过程中可以表达粘滞性，而在检波器波场逆时外推时则通过调整衰减系数的符号对正演过程的振幅损失进行补偿，使偏移成像更清晰. 本文将基于POD降维的有限元法［12-14］（以下简称POD有限元法）应用于粘滞声波方程的震源波场重构. 数值实验显示，相对于不使用震源波场重构策略的波场正演过程，POD有限元法能够以很少的计算时间代价减少存储量，且重构的震源波场具有较高精度，能很好地还原震源波场.

后文安排如下. 在第2节中我们给出粘滞声波方程的有限元法计算过程及其震源波场重构算法. 在第3节中我们给出基于POD有限元方法的震源波场重构的计算过程.在第4节中我们进行震源波场的数值模拟验证. 第5节为结论.

2 粘滞声波方程的有限元法

能够表达粘滞性的波动方程有多种. 从真实地震勘探成像出发，我们考虑如下简化的2维粘滞声波方程［11］的计算：

问题1 对2维粘滞声波方程的初边值问题，求u使得

c1ut（x，t）+c2utt（x，t）-Δu（x，t）=f（t）

x∈Ω，t∈（0，T），u（x，t）=0  x∈Ω，t∈［0，T］，

u（x，0）=φ0（x，y），ut（x，0）=φ1（x，y） x∈Ω（1）

其中ut=ut，utt=2ut2，Δu=2ux2+2uy2，c1=ψgc和c2=1c2为系数，其u=u（x，t）=u（x，y，t）为一个标量波场，x，y分别代表地表距离和深度距离，ft;xs为关于时间t的震源函数，xs为震源的位置，c代表波速，其中吸收系数g=πf0/cQ［15］，f0表示震源主频，Q是品质因子，ψ是常数，Ω瘙綆2为具有分片光滑边界Ω的区域，T是最后的计算时刻，源项为f（t），φ0（x，y）和φ1（x，y）为给定的初值函数. 为了进行波场模拟，以下我们令初值为0. 不同于一般的声波方程，如果ψ=1，方程（1）表示考虑介质衰减的正向延拓的震源波场. 如果ψ=-1，方程（1）表示考虑介质衰减补偿的反向延拓的检波器波场. 如果 ψ=0，方程（1）则退化为一般的声波方程. 因本文目的是研究粘滞波动方程的震源波场重构，以下的计算我们均令ψ=1.

为了由Galerkin原理导出问题（1）的变分形式，我们采用经典Sobolev空间［16］，即令U=H10（Ω）. 在方程两端同时乘任意测试函数v∈U，利用格林公式可得如下的关于问题（1）的变分形式.

问题2 对0

c1ut，v+c2utt，v+a（u，v）=（f，v），

v∈U，

u（x，t）=0， x∈Ω，t∈［0，T］，

u（x，0）=0， ut（x，0）=0， x∈Ω（2）

其中（·，·）是在L2（Ω）意义下的内积，a（u，v）=（SymbolQC@u，SymbolQC@v）.

令N为正整数，Δt=T/N为时间步长，进一步使用中心差分分别离散问题（2）中的二阶时间导数和一阶时间导数，即

untt=un+1-2un+un-1/Δt2，

unt=un+1-un-1/2Δt，

则问题（2）关于时间的半离散格式如下.

问题3 对Δt≤tn+1≤T-Δt，求un+1∈U，使得

c1un+1－un－1/2Δt，v+

c2un+1－2un+un－1/Δt2，v+

a（un，v）=（fn，v），v∈U，

u（x，0）=0， ut（x，0）=0， x∈Ω，

n=1，2，…，N－1（3）

在空间离散方法上，本文采用传统的有限元法进行离散. 令Ih是在Ω上的一致正则三角剖分［9］，有限元空间定义为

Uh={vh∈H10（Ω）∩C（Ω）：

vhK∈Pk（K），K∈Ih}，

其中Pk（K）是定义在单元K∈Ih上的小于等于k阶的多项式. 简单起见我们在接下来的计算中默认选取P1（K），即Lagrange线性元.

令unh∈Uh为问题3的解u在tn=nΔt（1≤n≤N-1）处的有限元近似解. 于是有如下关于有限元法的全离散格式.

问题4 求un+1h∈Uh，使得

c1un+1h－un－1h/2Δt，vh+

c2（un+1h－2unh+un－1h/

Δt2，vh）+a（unh，vh）=（fn，vh） vh∈Uh，u0h（x）=0， u0ht（x）=0，

n=1，2，…，N－1（4）

其中u0h（x）和u0ht（x）为时刻0处的初值，通过u0ht=u1h-u-1h/2Δt代入到问题（4）的离散格式中，消去u-1h就可以得到u1h的表达式.

接下来我们给出问题4的计算步骤. 如前所述，Nb维有限元子空间Uh可以近似Sobolev空间H10（Ω），其中Uh=spanφjNbj=1，Nb为剖分区域节点的数目. 于是，每一个解可以用该空间的一组基表示， 即用节点基函数的线性组合

un+1h（x，y）=∑Nbj=1un+1jφj（x，y）

表示有限元解，其中un+1j为待求时间层节点基函数的系数. 将其代入（4）式，选择vh=φi（x，y），i=1，…，Nb）为测试函数，整理可得如下的代数方程组：

c1M2Δt+c2MΔt2un+1j=2c2MΔt2-c2Aunj+

c1M2Δt-c2MΔt2un-1j+b→n，

n=1，2，…，N-1; j=1，2，…，Nb（5）

其中

A=aijNbi，j=1=∫ΩSymbolQC@φj·SymbolQC@φidxdyNbi，j=1，

M=mijNbi，j=1=∫ΩφjφidxdyNbi，j=1，

b→n=biNbi=1=∫ΩfφidxdyNbi=1

分别称为刚度矩阵、质量矩阵及载荷向量. 利用（5）式求得系数向量之后，再由un+1h（x，y）=∑Nbj=1un+1jφj（x，y）我们就得到有限元近似解.

下面我们给出基于有限元法的粘滞震源波场重构过程. 我们先利用（5）式计算出整个时间长度T的震源波场. 在震源波场的重构过程中，由于要逆向外推震源波场，我们只需在正演计算过程中存储所有的时间层边界值，以及最后两个时间层的波场值. 不同于式（4）和（5），用有限元法重构震源波场过程的数值格式是在时间上逆向求解问题（4），即，

问题5 求un-1h∈Uh，使得

c1un+1h－un－1h/2Δt，vh+ c2un+1h－2unh+un－1h/Δt2，vh+ a（unh，vh）=（fn，vh），vh∈Uh，uNhx，y=utmax（x，y）， uN－1hx，y=

utmax－Δt（x，y），n=N－1，…，1（6）

也就是说，我们要以最后两时间层的波场作为初值进行逆时外推，其中tmax为最大时间层，tmax-Δt为tmax时刻的上一时间层，utmax（x，y）为tmax时刻由格式（5）计算得到的波场值，umax-Δt为tmax-Δt时刻由格式（5）计算得到的波场值. 那么，相对于式（5）的逆过程就要求解线性方程组

c2MΔt2-c1M2Δtun-1j=2c2MΔt2-c2Aunj-

c1M2Δt+c2MΔt2un+1j+b→n，

n=N-1，…，1; j=1，2，…，Nb（7）

3 基于POD方法的粘滞震源波场重构

在本节中，我们将在式（5）和（7）基础上给出POD方法的计算过程，并将其应用于粘滞震源波场重构.  我们称该方法为震源波场重构的POD有限元法.

在实际逆时偏移成像中，随着炮点的增加将会重复计算式（7）许多次，当网格剖分细密、计算时间长时会产生巨大的计算量. 这应当归因于有限元方法的计算维度太大. 因此，为了让震源波场的重构计算更具效率，我们在震源波场重构中采用POD方法对有限元方法进行降维，以节省重构计算量，同时保持重构的震源波场较高的精度，提高逆时偏移成像效率.

POD方法是一种自由度较少且有足够精度的降维方法，其本质是在最小二乘意义下寻找已知数据的一组正交基（称为POD基），即求解一组已知数据的最优逼近. Luo等［9］将POD方与一些偏微分方程的数值解相结合，把高维模型降为低维模型，极大地减少了计算量. 此外，Luo等还给出了双曲方程基于POD降维的有限元FE方法［17］，Tan等［8］ 则利用POD方法对声波方程的震源波场进行了重构. 综上，基于POD有限元方法的粘声波方程震源波场重构算法的实现过程如下.

步骤1 正演模拟粘滞声波波场，即先利用式（5）计算波场，得到震源波场ujh（x，y）在时间点tj=jΔt，j=0，1，…，J的时间序列，然后以固定的采样间隔提取L列震源波场值，得到采样矩阵Us=uw1huw2h…uwLh， 其中wi=1，2，…，L是列序号，通常取L<5.

步骤2 构造瞬像矩阵As=AijL×L，其中Aij=SymbolQC@uih，SymbolQC@ujh1L，（·，·）为通常的L2内积.

步骤3 求解瞬像矩阵的特征值和特征向量，即As大于零的特征值λ1≥λ2≥…≥λl>0和对应的特征向量vj=aj1，aj2，…，ajLT， j=1，2，…，l.

步骤4 选取POD基的数目. 在实际计算过程中，可以使用ζ=∑di=1λi/∑li=1λi≥0.99来选取POD基的数目记为Np.

步骤5 计算POD基的表达式，计算公式为

ψj=∑Li=1ajiuih/Lλj（j=1，2，…，Np）.

步骤6 求解POD有限元格式.

下面我们给出粘滞震源波场重构的POD有限元格式. 前5个步骤构造了POD基，这样就可以由POD基张成一个POD空间

Ud=spanψ1，ψ2，…，ψNp.

根据POD方法的相关理论［9］

un+1d（x，y）=∑Npj=1pn+1jψj（x，y），

n=1，2，…，N-1

可以近似有限元解un+1h（x，y），其中un+1d（x，y）称为POD有限元格式的POD解，pn+1j为POD基函数的系数. 记由POD基组成的矩阵为POD矩阵，即 P=ψ1ψ2…ψNp.于是得到如下关于式（6）和（7）的POD有限元：

问题6 求un-1d∈Ud，使得

c1un+1d－un－1d/2Δt，vd+

c2un+1d－2und+un－1d/Δt2，vd +a（und，vd）=（fn，vd），vd∈Ud，uNdx，y=utmax（x，y）， uN－1dx，y=

umax－Δt（x，y）， n=N－1，…，1（8）

其中最后两个时间层的POD解就用实际存储的tmax时刻和tmax-Δt时刻的波场值umax（x，y）和umax-Δt（x，y）来近似.

将POD解的表达式代入（8），整理后待求的未知量为pn-1j，导出的线性方程组为

c2MdΔt2-c1Md2Δtpn-1j=2c2MdΔt2-c2Adpnj-

c1Md2Δt+c2MdΔt2pn+1j+bdn，  n=N-1，…，1; j=1，2，…，Np（9）

其中POD方法的刚度矩阵、质量矩阵以及载荷向量分别和原FE方法有如下关系：

Ad=aijNpi，j=1=∫ΩSymbolQC@ψj·SymbolQC@ψidxdyNpi，j=1=

PTAP，

Md=mijNpi，j=1=∫ΩψjψidxdyNpi，j=1=

PTMP，

bdn=biNpi=1=∫ΩfψidxdyNpi=1=PTb→n.

由（9）式计算得到pn-1j后，利用POD基的线性组合可得POD有限元解un-1d，n=N-1，…，1，进而逆推得到所有波场值. 可见，在使用有限元（6）和（7）对粘滞声波震源波场重构的过程中，每个时间层求解的未知量个数为Nb，且Nb往往会随着网格加密变得非常大. 另一方面，在使用POD有限元（式（8）和（9））对粘滞声波震源波场重构的过程中，每个时间层求解的未知量个数仅为Np，且Np

4 数值算例

在本节中，我们分别在粘滞声波震源波场重构中对有限元法和POD有限元法进行测试. 我们以2维均匀介质模型为例，计算区域为0≤x，y≤1 km，区域剖分方式如图1所示，剖分为131 072个三角单元，网格节点数为66 049，震源选为Ricker子波，主频f0=20Hz，函数形式为

f（t）=sin2πf0texp-π2f20t24.

我们将震源放置在计算区域的中心即（0.5 km，05 km）处，时间步长Δt=0.0005 s，最大记录时间为0.15 s，品质因子Q=30.0，波速c=4000 m/s.

我们分别记录时间为0.05、0.10和0.15 s处的波场快照，图2为使用有限元法正演计算震源波场得到的波场快照时间序列. 在正演计算中，我们每隔15个采样点存储震源波场，并由这些采样波场构造POD基和POD矩阵. 由POD有限元的定义可知，构成的瞬像矩阵是一个20阶的矩阵，即L=20.

图3为使用有限元法重构震源波场得到的波场快照时间序列，可见二者并没有明显区别. 由于重构震源波场是正演震源波场的逆过程，因而理论上也应当这样. 但是，由于存在计算误差，重构后的震源波场相对正演的震源波场有一定误差，不过非常小. 图4为使用POD有限元法重构的震源波场，与图1相比也并无明显差异.

为了得到POD有限元方法重构的震源波场相对于正演模拟波场的最大误差，即max（u-ud），我们令u为通过式（5）计算得到的解，ud为通过式（9）计算得到的解. 我们记录T=0.075 s时的波场快照并提取不同的地表位置0.31，0.42和0625 km处的波形图进行比较. 如图5所示，POD有限元法所重构的波形图和正演模拟的波形图十分吻合，重构波场的最大波场误差约为0.006 25，在整个T时间内POD有限元法重构震源波场的计算时间仅为9.25 s，而有限元法的计算时间却达到了242.36 s. 由此可见，POD有限元法重构后的震源波场在保持较高精度的同时节省计算时间，而且这个优势随着计算时间的增加变得更加明显. 另一方面，如果不考虑震源波场重构，基于有限元法的正演波场模拟需要的计算机存储为100%，采用有限元法重构震源波场后的存储量约为全存储的2.23%，而采用POD有限元法重构震源波场的存储量则约为全存储的8.21%. 表1对比了图3和图4所示的两种方法的最大波场误差、计算效率以及存储量.

5 结 论

在逆时偏移成像中，震源波场重构技术能够大幅降低逆时偏移成像对波场存储量的需求. 针对传统的震源波场重构方法将引入额外计算复杂度的问题，我们提出了POD有限元法并将其应用于粘滞声波方程的震源波场重构. 数值算例显示，与经典有限元法相比，该POD有限元法在保持较高精度的同时可以节省存储空间，大大减少重构震源波场的时间.

参考文献：

［1］ Zhang Y， Zhang P， Zhang H. Compensating for visco-acoustic effects in reverse-time migration ［M］.  Tulsa： Society of Exploration Geophysicists， 2010.

［2］ Fei T W， Luo Y， Yang J， et al. Removing false images in reverse time migration： the concept of de-primary ［J］. Geophysics， 2015， 80： 237.

［3］ Dussaud E， Symes W W， Williamson P， et al. Computational strategies for reverse-time migration ［M］. Tulsa： Society of Exploration Geophysicists， 2008.

［4］ 唐晨， 王德利. 基于源波场重建与波场分解的逆时偏移［J］. 世界地质， 2012， 31： 803.

［5］ 李松龄. 基于波场重建的时域全波形反演［J］. 中国锰业， 2018， 36： 8.

［6］ 王烈衡， 许学军. 有限元方法的数学基础［M］. 北京： 科学出版社， 2004.

［7］ Ciarlet P G. The finite element method for elliptic problems ［M］. Philadelphia： SIAM， 2002.

［8］ Tan W， Li R， Wu B， et al. Seismic source wavefield reconstruction based on proper orthogonal decomposition for finite element modeling ［C］// SEG 2021 Workshop： 4th International Workshop on Mathematical Geophysics： Traditional & Learning， Virtual， 17-19 December 2021. Tulsa： Society of Exploration Geophysicists， 2022.

［9］ Luo Z， Chen G. Proper orthogonal decomposition methods for partial differential equations ［M］. New York： Academic Press， 2018.

［10］ Yang Z， Liu Y， Ren Z. Comparisons of visco-acoustic wave equations ［J］. J Geophys Eng， 2014， 11： 025004.

［11］ Deng F，McMechan G A. True-amplitude prestack depth migration ［J］. Geophysics， 2007， 72： 155.

［12］ Luo Z， Li H， Zhou Y， et al. A reduced finite element formulation based on POD method for two-dimensional solute transport problems ［J］. J  Math Anal Appl， 2012， 385： 371.

［13］ Luo Z D， Chen J， Sun P， et al. Finite element formulation based on proper orthogonal decomposition for parabolic equations ［J］. Sci China Ser A： Math， 2009， 52： 585.

［14］ Luo Z， Zhou Y， Yang X. A reduced finite element formulation based on proper orthogonal decomposition for Burgers equation ［J］. Appl Numer Math， 2009， 59： 1933.

［15］ Futterman W I. Dispersive body waves ［J］. J Geophys Res， 1962， 67： 5279.

［16］ Adams R A， Fournier J J F. Sobolev spaces ［M］. Amsterdam： Elsevier， 2003.

［17］ LuoZ D， Ou Q L， Wu J R. A reduced FE formulation based on POD method for hyperbolic equations ［J］. Acta Math Sci， 2012， 32： 1997.

引用本文格式：

中 文：  宋淳， 王璐， 冯民富. 基于POD有限元法的粘滞声波方程震源波场重构［J］. 四川大学学报：  自然科学版， 2023， 60：  031005.

英 文：  Song C，Wang L，Feng M F. POD finite elementmethod for source wave field  reconstruction of viscous acoustic wave equations ［J］. J Sichuan Univ：  Nat Sci Ed， 2023， 60：  031005.