宋淳 王璐 馮民富

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

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

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

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

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

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

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

能夠表達粘滯性的波動方程有多種. 從真實地震勘探成像出發(fā)，我們考慮如下簡化的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為具有分片光滑邊界Ω的區(qū)域，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）

在空間離散方法上，本文采用傳統(tǒng)的有限元法進行離散. 令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為剖分區(qū)域節(jié)點的數目. 于是，每一個解可以用該空間的一組基表示， 即用節(jié)點基函數的線性組合

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

表示有限元解，其中un+1j為待求時間層節(jié)點基函數的系數. 將其代入（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方法的粘滯震源波場重構

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

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

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

步驟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 數值算例

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

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

我們將震源放置在計算區(qū)域的中心即（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為使用有限元法重構震源波場得到的波場快照時間序列，可見二者并沒有明顯區(qū)別. 由于重構震源波場是正演震源波場的逆過程，因而理論上也應當這樣. 但是，由于存在計算誤差，重構后的震源波場相對正演的震源波場有一定誤差，不過非常小. 圖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有限元法重構后的震源波場在保持較高精度的同時節(jié)省計算時間，而且這個優(yōu)勢隨著計算時間的增加變得更加明顯. 另一方面，如果不考慮震源波場重構，基于有限元法的正演波場模擬需要的計算機存儲為100%，采用有限元法重構震源波場后的存儲量約為全存儲的2.23%，而采用POD有限元法重構震源波場的存儲量則約為全存儲的8.21%. 表1對比了圖3和圖4所示的兩種方法的最大波場誤差、計算效率以及存儲量.

5 結 論

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

參考文獻：

［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］. 中國錳業(yè)， 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，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.