宋淳 王璐 馮民富

疊前逆時偏移是地震勘探中一種流行的地下結(jié)構(gòu)成像方法，其成像條件需要同時刻的震源波場值與檢波器波場值. 這在實(shí)際計(jì)算中就需要把正演模擬的所有時刻的震源波場數(shù)據(jù)全部存儲下來，存儲需求大. 雖然震源波場重構(gòu)技術(shù)可以降低對于波場數(shù)據(jù)的存儲需求，但會引入額外的計(jì)算復(fù)雜度. 為解決這個問題，本文提出了POD有限元法，并將其應(yīng)用于粘滯震源波場重構(gòu).這里的本征正交分解（Proper Orthogonal Decomposition， POD）方法是一種降維方法，能夠在降低數(shù)據(jù)量的同時提供足夠的計(jì)算精度.?? 數(shù)值算例顯示，該方法比傳統(tǒng)的有限元法更節(jié)省存儲空間，能夠加快重構(gòu)速度.

震源波場重構(gòu)； 粘滯聲波方程； 有限元法； 本征正交分解

O241.82A2023.031005

收稿日期： 2022-05-18

基金項(xiàng)目： 國家自然科學(xué)基金（11971337）

作者簡介： 宋淳 （2000-）， 男， 碩士研究生， 主要研究方向?yàn)橛?jì)算數(shù)學(xué). 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］是一種流行的復(fù)雜地質(zhì)體成像手段. 該方法可分為三部分，即震源波場的正向延拓，檢波器波場的逆時延拓以及成像條件. 其中，成像條件一般采用互相關(guān)成像，因而需要震源波場與檢波器波場同一時刻的波場值. 但是，由于震源波場外推和檢波器波場外推在時間上并不同步，往往需要將所有時間點(diǎn)的震源波場都存儲下來才能成像，進(jìn)而導(dǎo)致計(jì)算機(jī)存儲量劇增. 為解決這個問題，Dussaud等［3］首先提出了震源波場重構(gòu)法. 該方法先正向外推一次震源波場，然后再進(jìn)行一次震源波場的逆時外推，以保證震源波場和檢波器波場外推在時間上同步，減少存儲開銷.

經(jīng)過不斷地改進(jìn)，目前重構(gòu)震源波場方法已有多種［4］. 常見的重構(gòu)方法基于有限差分［5］. 有限差分法適合處理規(guī)則區(qū)域問題，且編程簡單，因而得到廣泛應(yīng)用，但它不能處理復(fù)雜地質(zhì)結(jié)構(gòu)及復(fù)雜邊界條件. 與有限差分法不同，有限元法［6，7］具有網(wǎng)格剖分靈活、對復(fù)雜結(jié)構(gòu)區(qū)域容易求解等優(yōu)勢，且基于變分原理時該方法還可以處理復(fù)雜邊界條件，在重構(gòu)震源波場時有潛在的優(yōu)勢. 其中，有限元的邊界值法是重構(gòu)震源波場的常用方法［8］. 該方法在震源波場正演過程中只需存儲最后兩個時間層的波場值作為震源波場重構(gòu)的初值. 然而，經(jīng)典的有限元法計(jì)算量大. 隨著網(wǎng)格剖分加密、計(jì)算時間增加，該方法需要求解大規(guī)模線性方程組，從而引入額外的計(jì)算復(fù)雜度，這個缺點(diǎn)就愈加明顯.

POD方法是一種能夠提供足夠精度而計(jì)算自由度較少的降維方法. 通過構(gòu)造POD基和低維POD空間，基于POD的數(shù)值方法的計(jì)算自由度遠(yuǎn)遠(yuǎn)小于有限元方法，因而能夠簡化計(jì)算、節(jié)省計(jì)算時間. 例如，Luo等［9］將POD方法推廣應(yīng)用于偏微分方程的數(shù)值計(jì)算格式，包括基于POD降維的有限差分法、有限元法、有限體積法等. 大量研究表明，將POD方法應(yīng)用于經(jīng)典有限元法會在保持有限元法優(yōu)勢的同時大大節(jié)省計(jì)算時間.

大多數(shù)疊前逆時偏移方法基于聲波方程，該方程模擬了理想情況下的地震波傳播規(guī)律. 然而，真實(shí)的地下介質(zhì)具有粘滯性，常規(guī)的聲波偏移方法并不能正確地處理地下結(jié)構(gòu)成像，因而研究粘滯聲波方程的數(shù)值解法是有意義的. 值得注意的是，能夠表達(dá)粘滯性的波動方程有多種［10］，其中Deng和McMechan［11］在聲波方程中添加波場的一階時間導(dǎo)數(shù)來調(diào)整振幅，使方程在波場正演過程中可以表達(dá)粘滯性，而在檢波器波場逆時外推時則通過調(diào)整衰減系數(shù)的符號對正演過程的振幅損失進(jìn)行補(bǔ)償，使偏移成像更清晰. 本文將基于POD降維的有限元法［12-14］（以下簡稱POD有限元法）應(yīng)用于粘滯聲波方程的震源波場重構(gòu). 數(shù)值實(shí)驗(yàn)顯示，相對于不使用震源波場重構(gòu)策略的波場正演過程，POD有限元法能夠以很少的計(jì)算時間代價減少存儲量，且重構(gòu)的震源波場具有較高精度，能很好地還原震源波場.

后文安排如下. 在第2節(jié)中我們給出粘滯聲波方程的有限元法計(jì)算過程及其震源波場重構(gòu)算法. 在第3節(jié)中我們給出基于POD有限元方法的震源波場重構(gòu)的計(jì)算過程.在第4節(jié)中我們進(jìn)行震源波場的數(shù)值模擬驗(yàn)證. 第5節(jié)為結(jié)論.

2 粘滯聲波方程的有限元法

能夠表達(dá)粘滯性的波動方程有多種. 從真實(shí)地震勘探成像出發(fā)，我們考慮如下簡化的2維粘滯聲波方程［11］的計(jì)算：

問題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為系數(shù)，其u=u（x，t）=u（x，y，t）為一個標(biāo)量波場，x，y分別代表地表距離和深度距離，ft;xs為關(guān)于時間t的震源函數(shù)，xs為震源的位置，c代表波速，其中吸收系數(shù)g=πf0/cQ［15］，f0表示震源主頻，Q是品質(zhì)因子，ψ是常數(shù)，Ω瘙綆2為具有分片光滑邊界Ω的區(qū)域，T是最后的計(jì)算時刻，源項(xiàng)為f（t），φ0（x，y）和φ1（x，y）為給定的初值函數(shù). 為了進(jìn)行波場模擬，以下我們令初值為0. 不同于一般的聲波方程，如果ψ=1，方程（1）表示考慮介質(zhì)衰減的正向延拓的震源波場. 如果ψ=-1，方程（1）表示考慮介質(zhì)衰減補(bǔ)償?shù)姆聪蜓油氐臋z波器波場. 如果 ψ=0，方程（1）則退化為一般的聲波方程. 因本文目的是研究粘滯波動方程的震源波場重構(gòu)，以下的計(jì)算我們均令ψ=1.

為了由Galerkin原理導(dǎo)出問題（1）的變分形式，我們采用經(jīng)典Sobolev空間［16］，即令U=H10（Ω）. 在方程兩端同時乘任意測試函數(shù)v∈U，利用格林公式可得如下的關(guān)于問題（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（Ω）意義下的內(nèi)積，a（u，v）=（SymbolQC@u，SymbolQC@v）.

令N為正整數(shù)，Δt=T/N為時間步長，進(jìn)一步使用中心差分分別離散問題（2）中的二階時間導(dǎo)數(shù)和一階時間導(dǎo)數(shù)，即

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

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

則問題（2）關(guān)于時間的半離散格式如下.

問題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）

在空間離散方法上，本文采用傳統(tǒng)的有限元法進(jìn)行離散. 令I(lǐng)h是在Ω上的一致正則三角剖分［9］，有限元空間定義為

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

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

其中Pk（K）是定義在單元K∈Ih上的小于等于k階的多項(xiàng)式. 簡單起見我們在接下來的計(jì)算中默認(rèn)選取P1（K），即Lagrange線性元.

令unh∈Uh為問題3的解u在tn=nΔt（1≤n≤N-1）處的有限元近似解. 于是有如下關(guān)于有限元法的全離散格式.

問題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的表達(dá)式.

接下來我們給出問題4的計(jì)算步驟. 如前所述，Nb維有限元子空間Uh可以近似Sobolev空間H10（Ω），其中Uh=spanφjNbj=1，Nb為剖分區(qū)域節(jié)點(diǎn)的數(shù)目. 于是，每一個解可以用該空間的一組基表示， 即用節(jié)點(diǎn)基函數(shù)的線性組合

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

表示有限元解，其中un+1j為待求時間層節(jié)點(diǎn)基函數(shù)的系數(shù). 將其代入（4）式，選擇vh=φi（x，y），i=1，…，Nb）為測試函數(shù)，整理可得如下的代數(shù)方程組：

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

分別稱為剛度矩陣、質(zhì)量矩陣及載荷向量. 利用（5）式求得系數(shù)向量之后，再由un+1h（x，y）=∑Nbj=1un+1jφj（x，y）我們就得到有限元近似解.

下面我們給出基于有限元法的粘滯震源波場重構(gòu)過程. 我們先利用（5）式計(jì)算出整個時間長度T的震源波場. 在震源波場的重構(gòu)過程中，由于要逆向外推震源波場，我們只需在正演計(jì)算過程中存儲所有的時間層邊界值，以及最后兩個時間層的波場值. 不同于式（4）和（5），用有限元法重構(gòu)震源波場過程的數(shù)值格式是在時間上逆向求解問題（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）

也就是說，我們要以最后兩時間層的波場作為初值進(jìn)行逆時外推，其中tmax為最大時間層，tmax-Δt為tmax時刻的上一時間層，utmax（x，y）為tmax時刻由格式（5）計(jì)算得到的波場值，umax-Δt為tmax-Δt時刻由格式（5）計(jì)算得到的波場值. 那么，相對于式（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方法的粘滯震源波場重構(gòu)

在本節(jié)中，我們將在式（5）和（7）基礎(chǔ)上給出POD方法的計(jì)算過程，并將其應(yīng)用于粘滯震源波場重構(gòu).? 我們稱該方法為震源波場重構(gòu)的POD有限元法.

在實(shí)際逆時偏移成像中，隨著炮點(diǎn)的增加將會重復(fù)計(jì)算式（7）許多次，當(dāng)網(wǎng)格剖分細(xì)密、計(jì)算時間長時會產(chǎn)生巨大的計(jì)算量. 這應(yīng)當(dāng)歸因于有限元方法的計(jì)算維度太大. 因此，為了讓震源波場的重構(gòu)計(jì)算更具效率，我們在震源波場重構(gòu)中采用POD方法對有限元方法進(jìn)行降維，以節(jié)省重構(gòu)計(jì)算量，同時保持重構(gòu)的震源波場較高的精度，提高逆時偏移成像效率.

POD方法是一種自由度較少且有足夠精度的降維方法，其本質(zhì)是在最小二乘意義下尋找已知數(shù)據(jù)的一組正交基（稱為POD基），即求解一組已知數(shù)據(jù)的最優(yōu)逼近. Luo等［9］將POD方與一些偏微分方程的數(shù)值解相結(jié)合，把高維模型降為低維模型，極大地減少了計(jì)算量. 此外，Luo等還給出了雙曲方程基于POD降維的有限元FE方法［17］，Tan等［8］ 則利用POD方法對聲波方程的震源波場進(jìn)行了重構(gòu). 綜上，基于POD有限元方法的粘聲波方程震源波場重構(gòu)算法的實(shí)現(xiàn)過程如下.

步驟1 正演模擬粘滯聲波波場，即先利用式（5）計(jì)算波場，得到震源波場ujh（x，y）在時間點(diǎn)tj=jΔt，j=0，1，…，J的時間序列，然后以固定的采樣間隔提取L列震源波場值，得到采樣矩陣Us=uw1huw2h…uwLh， 其中wi=1，2，…，L是列序號，通常取L<5.

步驟2 構(gòu)造瞬像矩陣As=AijL×L，其中Aij=SymbolQC@uih，SymbolQC@ujh1L，（·，·）為通常的L2內(nèi)積.

步驟3 求解瞬像矩陣的特征值和特征向量，即As大于零的特征值λ1≥λ2≥…≥λl>0和對應(yīng)的特征向量vj=aj1，aj2，…，ajLT， j=1，2，…，l.

步驟4 選取POD基的數(shù)目. 在實(shí)際計(jì)算過程中，可以使用ζ=∑di=1λi/∑li=1λi≥0.99來選取POD基的數(shù)目記為Np.

步驟5 計(jì)算POD基的表達(dá)式，計(jì)算公式為

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

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

下面我們給出粘滯震源波場重構(gòu)的POD有限元格式. 前5個步驟構(gòu)造了POD基，這樣就可以由POD基張成一個POD空間

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

根據(jù)POD方法的相關(guān)理論［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基函數(shù)的系數(shù). 記由POD基組成的矩陣為POD矩陣，即 P=ψ1ψ2…ψNp.于是得到如下關(guān)于式（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解就用實(shí)際存儲的tmax時刻和tmax-Δt時刻的波場值umax（x，y）和umax-Δt（x，y）來近似.

將POD解的表達(dá)式代入（8），整理后待求的未知量為pn-1j，導(dǎo)出的線性方程組為

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

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

其中POD方法的剛度矩陣、質(zhì)量矩陣以及載荷向量分別和原FE方法有如下關(guān)系：

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）式計(jì)算得到pn-1j后，利用POD基的線性組合可得POD有限元解un-1d，n=N-1，…，1，進(jìn)而逆推得到所有波場值. 可見，在使用有限元（6）和（7）對粘滯聲波震源波場重構(gòu)的過程中，每個時間層求解的未知量個數(shù)為Nb，且Nb往往會隨著網(wǎng)格加密變得非常大. 另一方面，在使用POD有限元（式（8）和（9））對粘滯聲波震源波場重構(gòu)的過程中，每個時間層求解的未知量個數(shù)僅為Np，且Np

4 數(shù)值算例

在本節(jié)中，我們分別在粘滯聲波震源波場重構(gòu)中對有限元法和POD有限元法進(jìn)行測試. 我們以2維均勻介質(zhì)模型為例，計(jì)算區(qū)域?yàn)?≤x，y≤1 km，區(qū)域剖分方式如圖1所示，剖分為131 072個三角單元，網(wǎng)格節(jié)點(diǎn)數(shù)為66 049，震源選為Ricker子波，主頻f0=20Hz，函數(shù)形式為

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

我們將震源放置在計(jì)算區(qū)域的中心即（0.5 km，05 km）處，時間步長Δt=0.0005 s，最大記錄時間為0.15 s，品質(zhì)因子Q=30.0，波速c=4000 m/s.

我們分別記錄時間為0.05、0.10和0.15 s處的波場快照，圖2為使用有限元法正演計(jì)算震源波場得到的波場快照時間序列. 在正演計(jì)算中，我們每隔15個采樣點(diǎn)存儲震源波場，并由這些采樣波場構(gòu)造POD基和POD矩陣. 由POD有限元的定義可知，構(gòu)成的瞬像矩陣是一個20階的矩陣，即L=20.

圖3為使用有限元法重構(gòu)震源波場得到的波場快照時間序列，可見二者并沒有明顯區(qū)別. 由于重構(gòu)震源波場是正演震源波場的逆過程，因而理論上也應(yīng)當(dāng)這樣. 但是，由于存在計(jì)算誤差，重構(gòu)后的震源波場相對正演的震源波場有一定誤差，不過非常小. 圖4為使用POD有限元法重構(gòu)的震源波場，與圖1相比也并無明顯差異.

為了得到POD有限元方法重構(gòu)的震源波場相對于正演模擬波場的最大誤差，即max（u-ud），我們令u為通過式（5）計(jì)算得到的解，ud為通過式（9）計(jì)算得到的解. 我們記錄T=0.075 s時的波場快照并提取不同的地表位置0.31，0.42和0625 km處的波形圖進(jìn)行比較. 如圖5所示，POD有限元法所重構(gòu)的波形圖和正演模擬的波形圖十分吻合，重構(gòu)波場的最大波場誤差約為0.006 25，在整個T時間內(nèi)POD有限元法重構(gòu)震源波場的計(jì)算時間僅為9.25 s，而有限元法的計(jì)算時間卻達(dá)到了242.36 s. 由此可見，POD有限元法重構(gòu)后的震源波場在保持較高精度的同時節(jié)省計(jì)算時間，而且這個優(yōu)勢隨著計(jì)算時間的增加變得更加明顯. 另一方面，如果不考慮震源波場重構(gòu)，基于有限元法的正演波場模擬需要的計(jì)算機(jī)存儲為100%，采用有限元法重構(gòu)震源波場后的存儲量約為全存儲的2.23%，而采用POD有限元法重構(gòu)震源波場的存儲量則約為全存儲的8.21%. 表1對比了圖3和圖4所示的兩種方法的最大波場誤差、計(jì)算效率以及存儲量.

5 結(jié) 論

在逆時偏移成像中，震源波場重構(gòu)技術(shù)能夠大幅降低逆時偏移成像對波場存儲量的需求. 針對傳統(tǒng)的震源波場重構(gòu)方法將引入額外計(jì)算復(fù)雜度的問題，我們提出了POD有限元法并將其應(yīng)用于粘滯聲波方程的震源波場重構(gòu). 數(shù)值算例顯示，與經(jīng)典有限元法相比，該P(yáng)OD有限元法在保持較高精度的同時可以節(jié)省存儲空間，大大減少重構(gòu)震源波場的時間.

參考文獻(xiàn)：

［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］. 世界地質(zhì)， 2012， 31： 803.

［5］ 李松齡. 基于波場重建的時域全波形反演［J］. 中國錳業(yè)， 2018， 36： 8.

［6］ 王烈衡， 許學(xué)軍. 有限元方法的數(shù)學(xué)基礎(chǔ)［M］. 北京： 科學(xué)出版社， 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有限元法的粘滯聲波方程震源波場重構(gòu)［J］. 四川大學(xué)學(xué)報：? 自然科學(xué)版， 2023， 60：? 031005.

英 文：? Song C，Wang L，F(xiàn)eng 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.