起伏地表無網(wǎng)格井中地震Q補(bǔ)償逆時(shí)偏移成像方法

摘要：井中地震震源能量弱，地層吸收衰減效應(yīng)的影響更強(qiáng)；復(fù)雜地表及復(fù)雜構(gòu)造對(duì)井周精細(xì)成像影響巨大，需要發(fā)展針對(duì)性的井中地震起伏地表衰減補(bǔ)償偏移成像方法。在徑向基函數(shù)有限差分方法原理基礎(chǔ)上，通過引入混合高斯立方基函數(shù)并設(shè)計(jì)起伏地表節(jié)點(diǎn)剖分策略，提出一種高精度的無網(wǎng)格徑向基函數(shù)起伏地表黏聲波數(shù)值模擬方法；將該方法應(yīng)用至改進(jìn)的基于Kelvin-Voigt模型的黏聲波方程，實(shí)現(xiàn)高效穩(wěn)定的井中地震起伏地表無網(wǎng)格Q補(bǔ)償逆時(shí)偏移成像。結(jié)果表明，提出的井中地震Q補(bǔ)償逆時(shí)偏移成像方法可以實(shí)現(xiàn)高精度的起伏地表井中地震衰減補(bǔ)償成像，有效恢復(fù)剖面能量、校正相位畸變，均衡剖面中的能量分布，提高成像質(zhì)量。

關(guān)鍵詞：井中地震； 起伏地表； 徑向基函數(shù)有限差分； Q補(bǔ)償逆時(shí)偏移

中圖分類號(hào)：P 631.4"" 文獻(xiàn)標(biāo)志碼：A

引用格式：魏國(guó)華，谷丙洛，段沛然，等.起伏地表無網(wǎng)格井中地震Q補(bǔ)償逆時(shí)偏移成像方法［J］.中國(guó)石油大學(xué)學(xué)報(bào)（自然科學(xué)版），2024，48（4）：68-79.

WEI Guohua， GU Bingluo， DUAN Peiran， et al. Q-compensated borehole seismic data reverse time migration with irregular topography based on mesh free method［J］. Journal of China University of Petroleum（Edition of Natural Science）， 2024，48（4）：68-79.

Q-compensated borehole seismic data reverse time migration with irregular topography based on mesh free method

WEI Guohua1，2， GU Bingluo1， DUAN Peiran1， ZHANG Shanshan1， LI Zhenchun1， KONG Qingfeng2

（1.School of Geosciences in China University of Petroleum （East China）， Qingdao 266580， China;

2.Shengli Geophysical Institute in SINOPEC， Dongying 257022， China）

Abstract：Borehole seismic theoretically offers advantages such as" high signal-to-noise ratio， wide frequency band， rich wavefield information， and strong reservoir identification capabilities， allowing for detailed reservoir imaging around the well. However， the source energy of borehole seismic is weak， resulting in stronger absorption attenuation effect in subsurface media compared to surface seismic. Furthermore， irregular topography and complex structure around the borehole significantly impact fine migration imaging. Therefore， it is essential to develop an attenuation compensation migration imaging method suited for irregular topography in borehole seismic. Using the principle of radius-basis-function finite-difference method， we establish a high-precision meshfree radius-basis-function finite-difference method with irregular topography of viscoacoustic wave equation. This method introduces a mixed Gaussian cubic basis-function and designs a node generation strategy for irregular topography. In addition， we apply this method to the improved viscoacoustic wave equation based on Kelvin-Voigt model to implement an efficient and stable Q-compensated reverse time migration. The results demonstrate that the proposed Q-compensated reverse time migration imaging method can achieve high-precision imaging of irregular topography for borehole seismic with attenuation compensation， effectively restoring the profile energy， correcting phase distortion， balancing energy distribution in the profile， and improving" imaging quality.

Keywords： borehole seismic; irregular topography; radius-basis-function finite-difference; Q-compensated reverse time migration

井中地震技術(shù)是指將檢波器布設(shè)于地表，震源位于數(shù)百米深的井內(nèi)激發(fā)的地震勘探方法。區(qū)別于常規(guī)勘探范圍廣、探測(cè)深度大的地面地震勘探，井中地震數(shù)據(jù)的采集、接收系統(tǒng)或能量激發(fā)系統(tǒng)靠近目標(biāo)地質(zhì)體，采集的數(shù)據(jù)包含豐富的目標(biāo)層信息，可針對(duì)特定油藏目標(biāo)進(jìn)行小范圍勘探，能克服復(fù)雜環(huán)境采集難題，減少震源施工困難，提高施工效率；井中地震波場(chǎng)信息豐富且干擾少，波場(chǎng)特征易于識(shí)別與處理，可實(shí)現(xiàn)反射波成像、直達(dá)波層析成像、轉(zhuǎn)換橫波成像等。鑒于以上優(yōu)勢(shì)，發(fā)展井中地震技術(shù)對(duì)于當(dāng)前的復(fù)雜油氣藏勘探具有重要的實(shí)用價(jià)值。近些年，油氣勘探靶區(qū)主要面臨劇烈起伏的地表及復(fù)雜的近地表結(jié)構(gòu)，常規(guī)基于靜校正處理的偏移成像結(jié)果發(fā)生畸變，成像精度難以滿足后續(xù)地質(zhì)地球物理解釋的需求。對(duì)此復(fù)雜問題，基于起伏地表的直接成像是最為有效的解決方法。目前，基于起伏地表的直接成像已經(jīng)發(fā)展出較多方法［1-7］，盡管這些方法應(yīng)用于實(shí)際生產(chǎn)并取得了一定的效果，但針對(duì)井中地震的起伏地表偏移成像方法依然有待深入研究。此外，實(shí)際地下介質(zhì)廣泛存在吸收衰減效應(yīng)，這對(duì)地震偏移成像也造成嚴(yán)重的影響，導(dǎo)致成像不準(zhǔn)，甚至失敗。因觀測(cè)系統(tǒng)特殊，介質(zhì)吸收衰減效應(yīng)對(duì)井中地震資料的影響要強(qiáng)于傳統(tǒng)地面地震資料，因此研發(fā)針對(duì)性的吸收衰減補(bǔ)償偏移成像方法對(duì)于井中地震技術(shù)的實(shí)際推廣應(yīng)用具有重要的理論及實(shí)際意義。Q補(bǔ)償逆時(shí)偏移（Q-RTM）基于雙程波動(dòng)方程，精度高，近年來逐漸成為業(yè)內(nèi)研究和發(fā)展的重點(diǎn)?；谒灭阅Ｐ筒煌l(fā)展出多種不同的Q-RTM成像方法［8-12］。盡管針對(duì)衰減補(bǔ)償逆時(shí)偏移成像方法的研究已經(jīng)較為深入，但針對(duì)井中地震的起伏地表衰減補(bǔ)償逆時(shí)偏移成像方法仍需進(jìn)一步深入研究。面對(duì)井中地震的井周構(gòu)造精細(xì)成像的實(shí)際需求，起伏地表黏聲逆時(shí)偏移的理論方法研究仍然不夠完善。筆者在常規(guī)徑向基函數(shù)有限差分方法原理基礎(chǔ)上，通過引入混合高斯立方基函數(shù)并設(shè)計(jì)起伏地表節(jié)點(diǎn)剖分策略，發(fā)展一種高精度的無網(wǎng)格徑向基函數(shù)起伏地表黏聲波數(shù)值模擬方法，并將所提出的數(shù)值模擬方法應(yīng)用至改進(jìn)的基于Kelvin-Voigt模型的黏聲波方程，以克服現(xiàn)有基于貼體網(wǎng)格的黏聲波方程數(shù)值求解方法對(duì)地形起伏的要求，實(shí)現(xiàn)高效穩(wěn)定的井中地震起伏地表無網(wǎng)格黏聲逆時(shí)偏移成像。

1 方法原理

1.1 Q補(bǔ)償逆時(shí)偏移基本原理

疊前逆時(shí)偏移主要包括3個(gè)步驟：

（1）震源波場(chǎng)正傳。具體可以表示為

LpS（x，t）=f（t）δ（x-xS）.（1）

式中，L為地震波場(chǎng)正傳算子；pS（x，t）為震源波場(chǎng)；f（t）表示震源時(shí)間函數(shù)，即震源子波；δ為脈沖函數(shù)；x為空間位置矢量；xS為震源空間位置矢量；t為地震波傳播時(shí)間。

（2）檢波點(diǎn)波場(chǎng)反傳。具體表示為

LTpR（x，t）=d（xR，t;xS）.（2）

式中，LT為地震波場(chǎng)延拓伴隨算子；pR（x，t）為檢波點(diǎn)波場(chǎng)；d為地震記錄；xR為檢波點(diǎn)空間位置矢量。

（3）應(yīng)用成像條件，獲得偏移成像剖面。基于常用的互相關(guān)成像條件，該步驟表示為

I（x）=∫tm0pS（x，t）pR（x，t）dt.（3）

式中，I（x）為偏移成像剖面；tm為地震記錄時(shí)長(zhǎng)。

圖1（a）為無衰減介質(zhì)地震波場(chǎng)延拓過程。由圖1可知，在無衰減介質(zhì)波場(chǎng)延拓過程中不存在吸收衰減效應(yīng)，因此直接應(yīng)用上述3個(gè)步驟即可實(shí)現(xiàn)對(duì)地下構(gòu)造的準(zhǔn)確成像。然而，對(duì)于衰減介質(zhì)，其地震波場(chǎng)延拓過程如圖1（b）所示，由圖可知，地震波在衰減介質(zhì)中傳播時(shí)存在明顯的吸收衰減效應(yīng)，地震波由震源經(jīng)反射面至檢波點(diǎn)處時(shí)，地震波場(chǎng)能量產(chǎn)生了顯著的衰減。經(jīng)過衰減的地震記錄dQ為

dQ=dexp（-αlD）exp（-αlU）.（4）

式中，lD和lU分別為下行波場(chǎng)和上行波場(chǎng)的傳播距離；α為衰減系數(shù)；exp（-αlD）為下行波在傳播路徑上的衰減；exp（-αlD）為上行波在傳播路徑上的衰減。

由圖1可知，對(duì)于衰減介質(zhì)，實(shí)現(xiàn)地下介質(zhì)構(gòu)造準(zhǔn)確成像需在波場(chǎng)延拓過程中對(duì)介質(zhì)吸收衰減效應(yīng)進(jìn)行補(bǔ)償。圖2（a）、（b）和（c）分別給出了無衰減介質(zhì)、衰減介質(zhì)及衰減介質(zhì)Q補(bǔ)償逆時(shí)偏移的原理示意圖（圖中，上標(biāo)A表示無衰減聲學(xué)介質(zhì)，V表示帶衰減黏聲介質(zhì)，C表示衰減補(bǔ)償）。地震波在衰減介質(zhì)中傳播時(shí)，由于介質(zhì)的吸收衰減作用，會(huì)導(dǎo)致地震波的能量衰減，傳播速度變慢，導(dǎo)致相位滯后。由圖2及式（4）可知，若直接對(duì)采集到的衰減地震數(shù)據(jù)實(shí)施逆時(shí)偏移，而不進(jìn)行衰減補(bǔ)償，其中的震源波場(chǎng)和檢波點(diǎn)波場(chǎng)均存在吸收衰減效應(yīng)，成像剖面將存在振幅失真，同相軸錯(cuò)位等問題。因此對(duì)于衰減介質(zhì)需在震源波場(chǎng)正向延拓和檢波點(diǎn)波場(chǎng)反向延拓過程中分別對(duì)地震波場(chǎng)的吸收衰減效應(yīng)進(jìn)行補(bǔ)償，也即Q補(bǔ)償逆時(shí)偏移。

基于Q補(bǔ)償逆時(shí)偏移基本流程，衰減補(bǔ)償后的震源波場(chǎng)

pCS（x，t）和檢波點(diǎn)波場(chǎng)

pCR（x，t）分別滿足：

pCS（x，t）=pS（x，t）exp（+αlD），（5）

pCR（x，t）=R（x，t）exp（+αl

U）.（6）

式中，R（x，t）表征衰減介質(zhì)中的上行波場(chǎng)。R（x，t）計(jì)算公式為

R（x，t）=pR（x，t）exp（-αlD）

exp（-αlU）.（7）

利用補(bǔ)償后的震源波場(chǎng)pCS（x，t）和檢波點(diǎn)波場(chǎng)pCR（x，t）可以實(shí)現(xiàn)衰減Q補(bǔ)償逆時(shí)偏移成像，具體為

IC（x）=∫tm0pCS（x，t）pCR（x，t）dt.（8）

將式（5）～（7）代入式（8）并整理可得

IC（x）=∫tm0pS（x，t）pR（x，t）dt.（9）

由式（9）可知，經(jīng)過補(bǔ)償后的偏移剖面IC（x）理論上等于無衰減介質(zhì)逆時(shí)偏移得到的偏移剖面I（x），即介質(zhì)的吸收衰減效應(yīng)得到完全補(bǔ)償。

1.2 無網(wǎng)格徑向基函數(shù)有限差分黏聲波方程數(shù)值模擬

基于Kelvin-Voigt彈性單元體模型，各向同性黏聲波方程可以表示為

1v22pt2=2p+ψω0Q2pt .（10）

式中，p=p（x，t）為黏聲波場(chǎng)；v（x）為介質(zhì)中聲波傳播速度；Q（x）為介質(zhì)品質(zhì)因子；2為拉普拉斯算子；ω0為參考角頻率;ψ為符號(hào)因子。式（10）中，ψ=1表示方程考慮介質(zhì)吸收衰減，ψ=-1表示方程考慮介質(zhì)吸收衰減補(bǔ)償。

式（10）右端第二項(xiàng)為振幅衰減項(xiàng)，但是對(duì)于介質(zhì)吸收衰減效應(yīng)中的相位頻散，該方程無法描述。結(jié)合標(biāo)準(zhǔn)線性體模型及常Q模型頻散關(guān)系的黏聲波動(dòng)方程可表示為

1v22pt2=1+1Qcos1γπγ22p+ψω0Q2pt .（11）

其中

γ≈arctan（1/Q）/π.

式中，γ為無量綱參數(shù)，通常滿足0lt;γlt;0.5。對(duì)比式（10）和（11）可知，式（11）不僅可以有效描述地震波場(chǎng)的振幅衰減，而且也可以描述其相位頻散。同時(shí)，式（11）僅僅包含常規(guī)的偏導(dǎo)數(shù)運(yùn)算，可利用時(shí)間-空間域有限差分進(jìn)行高效求解，計(jì)算效率較高。然而，常規(guī)有限差分法［13-15］難以適用于起伏地表。無網(wǎng)格法地震波數(shù)值模擬是近些年逐漸興起的一類數(shù)值模擬新方法，基于無網(wǎng)格的徑向基函數(shù)有限差分法能夠有效避免虛假反射，具有剖分節(jié)點(diǎn)隨模型參數(shù)變化、起伏界面處理靈活等優(yōu)勢(shì)［16］。本文中將基于無網(wǎng)格的徑向基函數(shù)有限差分法引入黏聲波動(dòng)方程數(shù)值模擬中，以便獲得高精度的黏聲地震波場(chǎng)。

無網(wǎng)格徑向基函數(shù)有限差分法是基于一組局部無序的離散節(jié)點(diǎn)構(gòu)造近似函數(shù)，實(shí)現(xiàn)偏微分方程求解，該問題本質(zhì)上屬于徑向基函數(shù)插值問題。該插值問題具體可以描述為：若給定Ω區(qū)域內(nèi)無序離散節(jié)點(diǎn)ΞΩ={x0，x1，…，xM}RD分布在函數(shù)F所在空間，則函數(shù)F可由徑向基函數(shù)φ插值表示，具體形式為

F（x）≈∑Nj=0ωjφ（x-xj）.（12）

式中，F(xiàn)（x）為插值函數(shù)；N為該插值問題所使用的徑向基函數(shù)的個(gè)數(shù)；φ（x-xj）為采用的徑向基函數(shù)；自變量x-xj又可表示為r，即r=x-xj，表征中心節(jié)點(diǎn)x與其鄰近節(jié)點(diǎn)xj之間的歐幾里得距離；ωj為第j個(gè)基函數(shù)對(duì)應(yīng)的權(quán)重系數(shù)。

基于式（12）所述的插值問題，根據(jù)無網(wǎng)格徑向基函數(shù)有限差分的基本原理，對(duì)于黏聲地震波場(chǎng)，任意中心點(diǎn)節(jié)點(diǎn)x0處的空間偏導(dǎo)數(shù)Ap可以表示為

Ap（x）x=x0=∑Mj=0ajp（xj）， ∑Mj=0aj=0.（13）

式中，A為空間偏導(dǎo)數(shù)算子；M為任意中心節(jié)點(diǎn)x0處空間差分計(jì)算模板中所需鄰近節(jié)點(diǎn)的個(gè)數(shù)；下標(biāo)j為參與中心點(diǎn)節(jié)點(diǎn)x0處差分運(yùn)算的鄰近差分節(jié)點(diǎn)號(hào)；xj為差分模板中第j個(gè)鄰近節(jié)點(diǎn)；aj為無網(wǎng)格徑向基函數(shù)有限差分系數(shù)。將式（13）代入式（12），差分系數(shù)aj可確定唯一解。

φ（x0-x0）φ（x1-x0）…φ（xM-x0）φ（x0-x1）φ（x1-x1）…φ（xM-x1）φ（x0-xM）φ（x1-xM）…φ（xM-xM）×

a0a1aM=

Aφ（x-x0）x=x0Aφ（x-x1）x=x0Aφ（x-xM）x=x0.（14）

表1中給出了5種常用的徑向基函數(shù)φ（r）［17-18］，其中ε為形變參數(shù)。研究表明，徑向基函數(shù)不同，相應(yīng)的無網(wǎng)格徑向基函數(shù)有限差分的精度及穩(wěn)定性也不同［19］。高斯基函數(shù)和逆多元二次函數(shù)在任意維空間都是正定的，因此其差分系數(shù)具有唯一解。多元二次函數(shù)不是正定函數(shù)，但是其廣義Fourier變換還是大于0，若滿足條件，其差分系數(shù)仍然收斂于唯一解。多諧樣條函數(shù)是分段平滑的，與其他基函數(shù)相比，能夠?yàn)榫€性系統(tǒng)提供相對(duì)好的條件數(shù)。Wendland函數(shù)則是一類最低次數(shù)的正定且緊支的徑向基函數(shù)。截至目前，地震勘探領(lǐng)域通常采用逆多元二次函數(shù)型基函數(shù)。

基于無網(wǎng)格徑向基函數(shù)有限差分法求解式（11），首先將式（11）中的時(shí)間二階偏導(dǎo)數(shù)進(jìn)行二階差分離散，可得

pn+1=2pn-pn-1+v2dt21+1Qcos1γπγ22p+v2dt2ψω0Q2.（15）

式中，上標(biāo)n表示當(dāng)前時(shí)刻；為黏聲波場(chǎng)p的時(shí)間一階偏導(dǎo)數(shù)。采用無網(wǎng)格徑向基函數(shù)有限差分算子對(duì)（15）式中空間導(dǎo)數(shù)進(jìn)行離散，具體為

2px2x=xi≈∑Mj=1aijpj，

2pz2x=xi≈∑Mj=1bijpj.（16）

式中，上標(biāo)和下標(biāo)i為當(dāng)前中心節(jié)點(diǎn)號(hào)；aij和bij分別為x方向和z方向二階空間偏導(dǎo)數(shù)在中心節(jié)點(diǎn)xi處的無網(wǎng)格徑向基函數(shù)有限差分系數(shù)。同理，對(duì)于采用式（16）進(jìn)行空間偏導(dǎo)數(shù)計(jì)算。基于式（16），式（15）可以表示為

pn+1=2pn-pn-1+v2dt21+1Qcos1γπγ2×

∑Mj=1aijpj+∑Mj=1bijpj+v2dt2ψω0Q∑Mj=1aijj+∑Mj=1bijj.（17）

盡管地震勘探領(lǐng)域通常采用逆多元二次函數(shù)型基函數(shù)實(shí)現(xiàn)無網(wǎng)格徑向基函數(shù)有限差分，但是其穩(wěn)定性強(qiáng)烈依賴模擬參數(shù)，且穩(wěn)定性較差。常見的徑向基函數(shù)，如高斯基函數(shù)，理論上提供一種頻譜收斂的函數(shù)空間逼近偏微分方程的數(shù)值解。當(dāng)形變參數(shù)ε趨于0時(shí)，高斯基函數(shù)收斂為多項(xiàng)式插值，此時(shí)高斯基函數(shù)近似平坦，在無網(wǎng)格近似中具有更高的精度。但在無網(wǎng)格徑向基函數(shù)有限差分方法中采用較小的形變參數(shù)會(huì)導(dǎo)致差分矩陣病態(tài)，因此亟需一種改進(jìn)的無網(wǎng)格徑向基函數(shù)有限差分方法在高精度的同時(shí)解決病態(tài)問題。Mishra等［20］提出一種混合高斯立方基函數(shù)，用于提高差分方程數(shù)值求解的穩(wěn)定性，并將此方法應(yīng)用于離散數(shù)據(jù)插值問題。該混合高斯立方基函數(shù)的基本思想是結(jié)合高斯函數(shù)和立方函數(shù)的優(yōu)點(diǎn)，在不增加計(jì)算量的同時(shí)增強(qiáng)差分系數(shù)矩陣的穩(wěn)定性?；旌细咚沽⒎交瘮?shù)的具體表達(dá)式為

φ（r）=αexp［-（εr）2］+βr3.（18）

式中，α和β為分別用于控制高斯基函數(shù)和立方基函數(shù)貢獻(xiàn)的權(quán)重。Mishra等進(jìn)一步將該函數(shù)應(yīng)用于徑向基函數(shù)偽譜法求解偏微分方程［21］，計(jì)算精度高，穩(wěn)定性好。為了保證無網(wǎng)格徑向基函數(shù)有限差分方法在黏聲波方程數(shù)值求解中的穩(wěn)定性和精度，本文中將混合高斯立方基函數(shù)引入黏聲波方程數(shù)值求解，以便解決無網(wǎng)格徑向基函數(shù)有限差分中的病態(tài)問題。

為了驗(yàn)證本文方法的穩(wěn)定性和有效性，采用速度為2000 m/s，Q為10的均勻介質(zhì)模型進(jìn)行測(cè)試，模型大小為1000 m×1000 m，節(jié)點(diǎn)間隔為10 m，時(shí)間步長(zhǎng)為1 ms，最大記錄時(shí)間為2.0 s。震源為20 Hz雷克子波，位于近（500 m， 20 m）處。圖3（a）給出了規(guī)則節(jié)點(diǎn)分布圖，圖3（b）為均勻的無網(wǎng)格節(jié)點(diǎn)分布圖，采用自適應(yīng)無網(wǎng)格節(jié)點(diǎn)剖分策略［20］生成，二者空間間隔均為10 m，其中圖3（a）進(jìn)行了抽稀顯示。相比之下，無網(wǎng)格節(jié)點(diǎn)剖分基于一組局部無序的離散節(jié)點(diǎn)，并無結(jié)構(gòu)網(wǎng)格約束，因此具有更高的靈活性。圖4為采用不同方法進(jìn)行聲波方程和黏聲方程模擬得到的點(diǎn)（500 m， 0 m）處的波形圖。由圖可見，黏聲地震記錄中振幅明顯衰減，同時(shí)波形之間存在明顯時(shí)移。

與常規(guī)基函數(shù)相比，采用混合高斯立方基函數(shù)進(jìn)行模擬的波形更加穩(wěn)定，有效抑制了常規(guī)基函數(shù)易產(chǎn)生的波形振蕩，在波前的尾部仍然穩(wěn)定傳播。由此可知，基于混合高斯立方基函數(shù)的徑向基函數(shù)有限差分法比常規(guī)的徑向基函數(shù)有限差分法具有更高的模擬精度，在黏聲方程模擬中具有更高的適應(yīng)性和穩(wěn)定性。

1.3 井中地震起伏地表無網(wǎng)格Q補(bǔ)償逆時(shí)偏移

無網(wǎng)格節(jié)點(diǎn)剖分是徑向基函數(shù)有限差分的關(guān)鍵步驟，也是實(shí)現(xiàn)無網(wǎng)格徑向基函數(shù)有限差分地震波數(shù)值模擬的前提。傳統(tǒng)有限差分法通常采用均一的網(wǎng)格間距，然而采用規(guī)則網(wǎng)格對(duì)復(fù)雜構(gòu)造以及含低速帶的模型進(jìn)行數(shù)值模擬，會(huì)導(dǎo)致局部波場(chǎng)采樣不足而造成明顯的數(shù)值頻散，或者導(dǎo)致局部過采樣影響計(jì)算效率。不同于傳統(tǒng)的規(guī)則網(wǎng)格、貼體網(wǎng)格或非結(jié)構(gòu)三角網(wǎng)格，無網(wǎng)格節(jié)點(diǎn)相互之間無需任何連接，可以在一系列非規(guī)則的離散節(jié)點(diǎn)上形成差分模板，進(jìn)而實(shí)現(xiàn)波動(dòng)方程求解。當(dāng)前，無網(wǎng)格節(jié)點(diǎn)剖分主流為基于靜電排斥理論的節(jié)點(diǎn)生成法，主要包括3種類型：迭代法、推進(jìn)法和快速節(jié)點(diǎn)生成法［22］。在迭代法中，通過最小能量算法［23］或Lloyd算法［24］對(duì)速度模型進(jìn)行節(jié)點(diǎn)剖分，基于排斥力理論，逐步調(diào)整節(jié)點(diǎn)分布。然而，此類方法在每次迭代中需較大的存儲(chǔ)量和計(jì)算量，節(jié)點(diǎn)生成效率低。在推進(jìn)法中，節(jié)點(diǎn)由模型邊界向某個(gè)方向推進(jìn)剖分，直到區(qū)域完全覆蓋，因此推進(jìn)法無需多次迭代。與迭代法相比，推進(jìn)法可以顯著提高節(jié)點(diǎn)生成的效率，且實(shí)現(xiàn)過程簡(jiǎn)單。盡管推進(jìn)法被廣泛采用，但此類方法難以處理復(fù)雜模型，實(shí)際應(yīng)用中有明顯的局限性。Fornberg等［22］提出了快速節(jié)點(diǎn)生成方法，旨在優(yōu)化推進(jìn)法的節(jié)點(diǎn)剖分過程，提高剖分效率。該方法已在徑向基函數(shù)有限差分中成功應(yīng)用［25-26］。Duan等［19］在此方法基礎(chǔ)上，利用頻散關(guān)系和穩(wěn)定性條件約束節(jié)點(diǎn)間隔，提出了一種自適應(yīng)節(jié)點(diǎn)剖分方法，實(shí)現(xiàn)了靈活穩(wěn)定的徑向基函數(shù)有限差分聲波數(shù)值模擬。

本文中基于自適應(yīng)節(jié)點(diǎn)剖分方法，通過預(yù)設(shè)起伏地表邊界節(jié)點(diǎn)，在合理節(jié)點(diǎn)半徑范圍內(nèi)設(shè)置內(nèi)部節(jié)點(diǎn)，利用節(jié)點(diǎn)控制方程調(diào)整節(jié)點(diǎn)疏密性，形成了起伏地表節(jié)點(diǎn)剖分方法，具體流程如下：

（1）采用自適應(yīng)節(jié)點(diǎn)剖分法，根據(jù)模型幾何邊界（包括起伏地表），獲得邊界節(jié)點(diǎn)坐標(biāo)和相應(yīng)速度模型參數(shù)；

（2）根據(jù)自適應(yīng)節(jié)點(diǎn)剖分法，由邊界節(jié)點(diǎn)向模型內(nèi)部進(jìn)行剖分，獲得模型內(nèi)部節(jié)點(diǎn)的初始坐標(biāo)（ξ，η）和相應(yīng)速度模型參數(shù)；

（3）根據(jù)自適應(yīng)節(jié)點(diǎn)剖分法確定節(jié)點(diǎn)最小半徑，刪除內(nèi)部節(jié)點(diǎn)間距小于節(jié)點(diǎn)最小半徑一半的節(jié)點(diǎn)；

（4）根據(jù)自適應(yīng)節(jié)點(diǎn)控制方程，將模型內(nèi)部節(jié)點(diǎn)從初始坐標(biāo)調(diào)整到新的坐標(biāo)（x，z）（二者可視為映射關(guān)系），一般需要迭代4次以上才獲得滿足精度要求的節(jié)點(diǎn)分布，得到各節(jié)點(diǎn)的坐標(biāo)和速度模型參數(shù)。

在步驟（4）中，引入自適應(yīng)網(wǎng)格理論［27］，在獲得模型內(nèi)部節(jié)點(diǎn)分布的基礎(chǔ)上進(jìn)一步調(diào)整節(jié)點(diǎn)分布情況，其控制方程如下：

b1xξξ+b2xξη+b3xηη+a1zξξ+a2zξη+a3zηη=-J22wwx ，

a1xξξ+a2xξη+a3xηη+c1zξξ+c2zξη+c3zηη=-J22wwz .（19）

其中

a1=-xηzη， a2=xξzη+xηzξ， a3=-xξzξ，

b1=z2η， b2=-2zξzη， b3=z2ξ，

c1=x2η， c2=-2xξxη， c3=x2ξ，

J=xξzη-xηzξ.

式中，xξξ、xξη、xηη、zξξ、zξη和zηη分別為二階導(dǎo)數(shù)2x/ξ2、2x/ξη、2x/η2、2z/ξ2、2z/ξη和2z/η2的簡(jiǎn)寫；w（x，z）=1+αv2x+v2z為權(quán)重函數(shù)，其中α為權(quán)重函數(shù)的參數(shù)，vx=v/x和vz=v/z表示模型速度在x和z方向的一階導(dǎo)數(shù)；J為Jacobi多項(xiàng)式；xξ、xη、zξ和zη為一階偏導(dǎo)數(shù)x/ξ、x/η、z/ξ和z/η的簡(jiǎn)寫。

通過中心差分展開式（19），利用高斯-賽德爾迭代法求解。

基于起伏地表剖分的節(jié)點(diǎn)，結(jié)合前文所述無網(wǎng)格徑向基函數(shù)有限差分方法，即可實(shí)現(xiàn)井中地震起伏地表無網(wǎng)格Q補(bǔ)償逆時(shí)偏移。對(duì)于井中地震，地震波在衰減介質(zhì)中的傳播過程與地面地震類似，同樣會(huì)受吸收衰減效應(yīng)影響，因此基于圖2（c）所示基本原理，井中地震起伏地表無網(wǎng)格Q補(bǔ)償逆時(shí)偏移的實(shí)現(xiàn)流程如下：

（1）獲取偏移所需的模型參數(shù)，包括速度和Q，準(zhǔn)備待偏移的共炮點(diǎn)道集。

（2）基于上述起伏地表節(jié)點(diǎn)剖分方法，利用步驟（1）中模型參數(shù)實(shí)施節(jié)點(diǎn)剖分，獲得計(jì)算域內(nèi)節(jié)點(diǎn)位置。

（3）基于步驟（2）所得節(jié)點(diǎn)信息，利用無網(wǎng)格徑向基函數(shù)有限差分方法正向求解帶有衰減補(bǔ)償?shù)酿ぢ暡ǚ匠?，獲得補(bǔ)償后的震源波場(chǎng)，此處帶有衰減補(bǔ)償?shù)酿ぢ暡ǚ匠虨?/p>

1v22pt2=1+1Qcos1γπγ22p-1ω0Q2pt+f.

式中，f=f（x，t）為震源函數(shù)。

（4）同樣基于步驟（2）所得節(jié)點(diǎn)信息，利用無網(wǎng)格徑向基函數(shù)有限差分方法反向求解帶有衰減補(bǔ)償?shù)酿ぢ暡ǚ匠?，獲得補(bǔ)償后的檢波點(diǎn)波場(chǎng)，此處帶有衰減補(bǔ)償?shù)酿ぢ暡ǚ匠虨?/p>

1v22pt2=1+1Qcos1γπγ22p-1ω0Q2pt+dQ.

（5）將各離散節(jié)點(diǎn)對(duì)應(yīng)的震源波場(chǎng)和檢波點(diǎn)波場(chǎng)進(jìn)行最小曲率插值，獲得模型空間域規(guī)則網(wǎng)格點(diǎn)上對(duì)應(yīng)的波場(chǎng)值。

（6）對(duì)所得的規(guī)則網(wǎng)格點(diǎn)上對(duì)應(yīng)的震源波場(chǎng)和檢波點(diǎn)應(yīng)用成像條件（式（8）），獲得成像剖面。

在Q-RTM成像中，震源波場(chǎng)和檢波點(diǎn)波場(chǎng)均涉及衰減補(bǔ)償，在補(bǔ)償過程中除了有效信號(hào)外，波場(chǎng)中的數(shù)值高波數(shù)噪聲同樣也會(huì)被放大，導(dǎo)致波場(chǎng)延拓不穩(wěn)定。本文中通過低通濾波器壓制數(shù)值高波數(shù)噪聲，保障波場(chǎng)延拓的穩(wěn)定性。

2 數(shù)值算例

2.1 BP氣云模型

為了驗(yàn)證本文方法在復(fù)雜強(qiáng)衰減介質(zhì)中應(yīng)用的有效性，本文中首先采用大小為9 km×5 km的國(guó)際標(biāo)準(zhǔn)BP氣云模型進(jìn)行偏移試驗(yàn)，圖5為該模型的速度及Q值，模型中節(jié)點(diǎn)半徑為10 m。該模型中間位置存在氣云，具有強(qiáng)烈的吸收衰減作用。模擬中邊界以及起伏地表均采用PML吸收邊界條件?？紤]到井中地震有效成像范圍有限，此處共設(shè)計(jì)8口垂直井，起始井位于1 km處，井間距為1 km。每口井布設(shè)41個(gè)震源，起始震源位于地表之下200 m，震源垂向間距為20 m。井位的布設(shè)可見圖5中的紅色線段。檢波點(diǎn)沿起伏地表布設(shè)，單邊最大偏移距為2.5 km，道間距為10 m。記錄時(shí)間長(zhǎng)度4 s，時(shí)間步長(zhǎng)0.5 ms，采用主頻為25 Hz的雷克子波作為震源時(shí)間函數(shù)。

從圖5中的速度模型可以看出，該模型的速度變化較大。

圖6為BP氣云模型的自適應(yīng)剖分節(jié)點(diǎn)分布。在節(jié)點(diǎn)分布圖中，節(jié)點(diǎn)間隔隨著速度變化，淺層速度較小時(shí)，節(jié)點(diǎn)間距相對(duì)較小，節(jié)點(diǎn)密度較大，隨著速度增加，節(jié)點(diǎn)間距也隨之逐步增大，節(jié)點(diǎn)密度減小。使用自適應(yīng)節(jié)點(diǎn)剖分技術(shù)進(jìn)行網(wǎng)格剖分，可以大幅減少有效差分所需的網(wǎng)格點(diǎn)數(shù)，提供更加合理的節(jié)點(diǎn)分布和更加平滑的速度過渡，避免速度界面兩側(cè)出現(xiàn)節(jié)點(diǎn)缺失或異常分布，同時(shí)有效避免節(jié)點(diǎn)的缺失或重疊，使節(jié)點(diǎn)分布更加均勻?；诖穗x散節(jié)點(diǎn)實(shí)施波場(chǎng)延拓，可以有效提升計(jì)算效率。

圖7為BP氣云模型的偏移剖面，其中圖7（a）為黏聲數(shù)據(jù)RTM剖面，也即未補(bǔ)償RTM剖面，圖7（b）為黏聲數(shù)據(jù)Q-RTM剖面，圖7（c）為聲波數(shù)據(jù)RTM剖面，也即參考剖面。由圖7對(duì)比可知，在強(qiáng)衰減區(qū)域及其下方，常規(guī)的聲波RTM不能夠補(bǔ)償介質(zhì)的吸收衰減效應(yīng)，由于衰減造成的相位滯后影響了地震波的走時(shí)，這種相位畸變沒有得到校正，導(dǎo)致成像剖面中的同相軸不能完全收斂，同相軸的成像紊亂，同時(shí)由于沒有進(jìn)行衰減補(bǔ)償，剖面振幅失真，同相軸的能量較弱且不均衡，成像剖面與參考剖面相比質(zhì)量較低。相比之下，Q-RTM可以有效補(bǔ)償介質(zhì)的吸收衰減效應(yīng)對(duì)偏移的影響，對(duì)衰減區(qū)域下方的同相軸進(jìn)行能量恢復(fù)和相位校正，提高成像質(zhì)量及分辨率，其剖面質(zhì)量與參考剖面相當(dāng)，其結(jié)果驗(yàn)證了本文方法在強(qiáng)衰減介質(zhì)中進(jìn)行井中成像的有效性。圖8展示了氣云下方構(gòu)造的局部成像結(jié)果，從圖中可以看出，在經(jīng)過衰減補(bǔ)償之后，氣云下方構(gòu)造的成像結(jié)果得到顯著改善，其能量得到恢復(fù)，相位得到校正，說明本文方法可以有效進(jìn)行衰減補(bǔ)償。

2.2 勝利墾71薄層互模型

為了進(jìn)一步驗(yàn)證本文方法在復(fù)雜構(gòu)造井中成像的適用性，基于勝利油田墾71研究區(qū)實(shí)際地質(zhì)狀況，建立了典型的薄互層模型，如圖9所示。該模型地表起伏，右側(cè)發(fā)育正斷層，地層之間夾雜細(xì)碎的薄互層，且含有油氣，具有較強(qiáng)的吸收衰減特征。該模型大小為6000 m×2500 m，所用空間網(wǎng)格大小為5 m×5 m。該模型共設(shè)計(jì)6口垂直井，起始井位于橫向1000 m處，井間距為800 m。每口井布設(shè)41個(gè)震源，起始震源位于地表，震源垂向間距為20 m。井位的布設(shè)可見圖9中的紅色線段。檢波點(diǎn)沿起伏地表布設(shè)，單邊最大偏移距為3.0 km，道間距為10 m。地震記錄時(shí)間長(zhǎng)度1.6 s，時(shí)間步長(zhǎng)0.3 ms，采用主頻為60 Hz的雷克子波作為震源時(shí)間函數(shù)。

圖10為勝利墾71薄互層模型的自適應(yīng)剖分節(jié)點(diǎn)分布。本文中提出的自適應(yīng)節(jié)點(diǎn)剖分技術(shù)可以根據(jù)速度大小自主確定節(jié)點(diǎn)間隔，在合理節(jié)點(diǎn)半徑范圍內(nèi)設(shè)置內(nèi)部節(jié)點(diǎn)，自動(dòng)調(diào)整節(jié)點(diǎn)疏密，由自適應(yīng)剖分節(jié)點(diǎn)分布圖也可以看出，深層高速區(qū)域的節(jié)點(diǎn)間距大于淺層的低速區(qū)域，因此深層高速區(qū)域的節(jié)點(diǎn)數(shù)量也少于淺層的低速區(qū)域，從而大大減少節(jié)點(diǎn)數(shù)量，有效提高波場(chǎng)延拓的計(jì)算效率。

圖11為基于勝利墾71薄互層模型數(shù)值模擬得到的單炮地震記錄，震源位于（3400 m，420 m）處，其中圖11（a）為聲波地震記錄，圖11（b）為黏聲波地震記錄。由圖11對(duì)比可知，與聲波地震記錄相比，黏聲波地震記錄中的波場(chǎng)由于地下介質(zhì)的吸收衰減，導(dǎo)致同相軸的能量顯著減弱。為了更細(xì)致地對(duì)比該模型在不同介質(zhì)中的地震波場(chǎng)，在圖12中展示了從圖11中抽取的第300道的單道波形對(duì)比，其中虛線為聲波單道地震記錄，實(shí)線為衰減介質(zhì)中的單道記錄。從圖12對(duì)比可以看出，兩者的波形相似，但與無衰減介質(zhì)中的波形相比，衰減介質(zhì)中的地震波能量發(fā)生了顯著的衰減，相位出現(xiàn)明顯滯后，說明地下介質(zhì)的吸收衰減會(huì)同時(shí)影響地震波的振幅能量和相位信息。該波形對(duì)比圖更加直觀地展示了黏聲介質(zhì)對(duì)地震波的影響，同時(shí)證明本文方法可以有效地模擬介質(zhì)吸收衰減效應(yīng)對(duì)地震波場(chǎng)的影響。

圖13為基于勝利墾71薄互層模型所得的0.45 s波場(chǎng)快照，震源位于（3400 m，420 m）處。由圖可知，在衰減的波場(chǎng)快照中，地震波能量沒有得到補(bǔ)償，同相軸能量微弱，波場(chǎng)能量不均衡，波前面的位置不準(zhǔn)確；在經(jīng)過衰減補(bǔ)償之后，地震波能量得到恢復(fù)，同相軸能量增強(qiáng)，波前面位置準(zhǔn)確，波場(chǎng)能量分布均衡。

圖14為勝利墾71薄互層模型的逆時(shí)偏移結(jié)果，其中圖14（a）是黏聲數(shù)據(jù)RTM剖面，即未補(bǔ)償?shù)腞TM剖面，圖14（b）是黏聲數(shù)據(jù)Q-RTM剖面，即衰減補(bǔ)償?shù)腞TM剖面，圖14（c）為聲波數(shù)據(jù)RTM剖面，也即參考剖面。對(duì)比圖14可知，與參考剖面相比，未補(bǔ)償?shù)腞TM剖面中能量沒有得到恢復(fù)，同相軸能量微弱，深層不能有效、清晰地成像，剖面整體的能量分布不均衡，與參考剖面的差異明顯，成像質(zhì)量較低。經(jīng)過衰減補(bǔ)償之后圖14（b）中的成像剖面能量得到有效補(bǔ)償和恢復(fù)，同相軸能量和聚焦性增強(qiáng)，剖面整體的能量分布均衡，深層反射界面的成像清晰，反射界面的成像位置準(zhǔn)確，同相軸的連續(xù)性增強(qiáng)，薄互層的構(gòu)造刻畫清晰，分辨率提高，與參考剖面的成像質(zhì)量接近，成像質(zhì)量較高。為了更加細(xì)致地體現(xiàn)補(bǔ)償RTM的成像效果，圖15中展示了深層薄互層區(qū)域的局部放大。從圖中可以看出：未補(bǔ)償?shù)某上窠Y(jié)果中，薄互層區(qū)域同相軸能量微弱，同相軸連續(xù)性不強(qiáng)，無法有效成像；而在經(jīng)過衰減補(bǔ)償?shù)某上窠Y(jié)果中，薄互層區(qū)域的同相軸能量得到恢復(fù)，其聚焦性和連續(xù)性均得到提升，成像清晰，分辨率高，與參考剖面幾乎一致。

圖16為勝利墾71薄互層模型逆時(shí)偏移結(jié)果的波數(shù)譜（kx、kz分別為波數(shù)矢量在x和z方向的分量）

，其中圖16（a）是黏聲數(shù)據(jù)RTM剖面的波數(shù)譜，圖16（b）是黏聲數(shù)據(jù)Q-RTM剖面的波數(shù)譜，圖16（c）為聲波數(shù)據(jù)RTM剖面的波數(shù)譜。對(duì)比圖16可知：未補(bǔ)償剖面的波數(shù)譜的波數(shù)成分較少，低波數(shù)及高波數(shù)成分都不夠豐富，與參考剖面波數(shù)譜相差較大，成像剖面分辨率較低；而經(jīng)過衰減補(bǔ)償之后，成像剖面中的波數(shù)成分更加豐富，低波數(shù)及高波數(shù)成分均得到補(bǔ)充，波數(shù)譜的帶寬范圍增大，與參考剖面波數(shù)譜相近，說明衰減補(bǔ)償可以有效提高成像剖面的分辨率。

3 結(jié) 論

（1）起伏地表自適應(yīng)節(jié)點(diǎn)剖分技術(shù)克服了現(xiàn)有基于貼體網(wǎng)格的數(shù)值求解方法對(duì)地形起伏的要求，可以完全適應(yīng)地表地形，同時(shí)根據(jù)速度自主確定節(jié)點(diǎn)間隔，自動(dòng)調(diào)整節(jié)點(diǎn)疏密，在高速區(qū)域節(jié)點(diǎn)間距大而低速區(qū)域節(jié)點(diǎn)間距小，提供更加合理的節(jié)點(diǎn)分布和更加平滑的過渡，使節(jié)點(diǎn)分布更加均勻，極大減少網(wǎng)格點(diǎn)數(shù)，有效提升計(jì)算效率。

（2）在改進(jìn)的Kelvin-Voigt模型的黏聲波方程基礎(chǔ)上，結(jié)合自適應(yīng)節(jié)點(diǎn)剖分技術(shù)，推導(dǎo)并實(shí)現(xiàn)了無網(wǎng)格井中地震起伏地表黏聲逆時(shí)偏移成像。模型試算表明了該成像方法可以實(shí)現(xiàn)高精度的起伏地表井中地震衰減補(bǔ)償成像，有效恢復(fù)剖面能量、校正相位畸變，均衡剖面中的能量分布，提高成像質(zhì)量。

（3）應(yīng)用自適應(yīng)節(jié)點(diǎn)剖分技術(shù)，實(shí)現(xiàn)逆時(shí)偏移成像的思想在未來處理起伏地表復(fù)雜構(gòu)造、井中地震、衰減補(bǔ)償成像等領(lǐng)域具有廣闊的應(yīng)用前景。

參考文獻(xiàn)：

［1］ JAGER C， HERTWECK T， SPINER M. True-amplitude Kirchh off migration from topography ［C/OL］. 73rd Annual International Meeting， SEG， Dallas， Texas， United States of America， October 26-29， 2003［2022-08-10］.https：//library.seg.org/doi/epdfplus/10.1190/1.1818089.

［2］ GRAY S H. Gaussian beam migration of common-shot records［J］. Geophysics， 2005， 70（4）： S71-S77.

［3］ 黃建平，袁茂林，李振春，等.雙復(fù)雜條件下非傾斜疊加精確束偏移方法及應(yīng)用Ⅰ：聲波方程［J］.地球物理學(xué)報(bào)，2015，58（1）：257-276.

HUANG Jianping， YUAN Maolin， LI Zhenchun， et al. The accurate beam migration method without slant stack under dual-complexity conditions and its application（I）： acoustic equation［J］. Chinese Journal of Geophysics， 2015，58（1）：257-276.

［4］ 程玖兵，馬在田，陶正喜，等.山前帶復(fù)雜構(gòu)造成像方法研究［J］.石油地球物理勘探，2006，41（5）：525-529.

CHENG Jiubing， MA Zaitian， TAO Zhengxi， et al. Imaging study of piedmont complex structures［J］. Oil Geophysical Prospecting， 2006，41（5）：525-529.

［5］ 劉紅偉，劉洪，李博，等.起伏地表疊前逆時(shí)偏移理論及GPU加速技術(shù)［J］.地球物理學(xué)報(bào)，2011，54（7）：1883-1892.

LIU Hongwei， LIU Hong， LI Bo， et al. Pre-stack reverse time migration for furred topography and GPU acceleration technology ［J］. Chinese Journal of Geophysics， 2011，54（7）：1883-1892.

［6］ 黃韻博，黃建平，李振春，等.基于循環(huán)對(duì)抗神經(jīng)網(wǎng)絡(luò)的快速最小二乘逆時(shí)偏移成像方法［J］.中國(guó)石油大學(xué)學(xué)報(bào)（自然科學(xué)版），2023，47（3）：55-61.

HUANG Yunbo， HUANG Jianping， LI Zhenchun， et al. Fast least-squares reverse time migration based on cycle-consistent generative adversarial network［J］. Journal of China University of Petroleum （Edition of Natural Science）， 2023，47（3）：55-61.

［7］ 慕鑫茹，黃建平，黎國(guó)龍，等.黏聲TTI介質(zhì)中純qP波逆時(shí)偏移成像方法［J］.中國(guó)石油大學(xué)學(xué)報(bào)（自然科學(xué)版），2023，47（2）：44-52.

MU Xinru， HUANG Jianping， LI Guolong， et al. Reverse time migration in viscoacoustic TTI media based on pure qP wave equation［J］. Journal of China University of Petroleum（Edition of Natural Science）， 2023，47（2）：44-52.

［8］ DENG F， MCMECHAN G A. Viscoelastic true-amplitude prestack reverse-time depth migration［J］. Geophysics， 2008，73（4）：S143-S155.

［9］ ZHU T Y， HU W Y， NING J Y. An efficient local operator-based Q-compensated reverse time migration algorithm with multistage optimization ［J］. Geophysics， 2018，83（3）：S249-S259.

［10］ ZHU T Y， HARRIS J M， BIONDI B.Q-compensated reverse-time migration ［J］. Geophysics， 2014，79（3）：S77-S87.

［11］ YAN H， LIU Y.Visco-acoustic prestack reverse-time migration based on the time-space domain adaptive high-order finite-difference method［J］. Geophysical Prospecting， 2013，61（5）：941-954.

［12］ ZHAO Y， MAO N， REN ZM. A stable and efficient approach of Q reverse time migration［J］. Geophysics， 2018，83（6）：S557-S567.

［13］ KELLY K R， WARD R W， TREITEL S， et al. Synthetic seismograms： a finite-difference approach［J］.Geophysics， 1976，41（1）：2-27.

［14］ MARFURT K J. Accuracy of finite-difference and finite-element modeling of the scalar and elastic wave equations［J］. Geophysics，1984，49（5）：533-549.

［15］ VIRIEUX J. P-SV wave propagation in heterogeneous media： velocity-stress finite-difference method［J］. Geophysics， 1986，51（4）：889-901.

［16］ 劉立彬，段沛然，張?jiān)沏y，等.基于無網(wǎng)格的地震波場(chǎng)數(shù)值模擬方法綜述［J］.地球物理學(xué)進(jìn)展，2020，35（5）：1815-1825.

LIU Libin， DUAN Peiran， ZHANG Yunyin， et al. Overview of mesh-free method of seismic forward numerical simulation［J］. Progress in Exploration Geophysics， 2000，35（5）：1815-1825.

［17］ BUHMANN M D. Radial basis functions： theory and implementations［M］. Cambridge： Cambridge University Press， 2003：99-136.

［18］ SU L. A radial basis function （RBF）-finite difference （FD） method for the backward heat conduction problem［J］. Applied Mathematics and Computation， 2019，354：232-247.

［19］ DUAN P R， GU B L， LI Z C， et al. An adaptive node-distribution method for radial-basis-function finite difference modeling with optimal shape parameter［J］. Geophysics， 2021，86（1）：T1-T18.

［20］ MISHRA P K， NATH S K， KOSEC G， et al. An improved radial basis-pseudospectral method with hybrid Gaussian-cubic kernels［J］. Engineering Analysis with Boundary Elements， 2017，80：162-171.

［21］ MISHRA P K， NATH S K， SEN M K， et al. Hybrid Gaussian-cubic radial basis functions for scattered data interpolation［J］. Computational Geosciences， 2018，22（5）：1203-1218.

［22］ FORNBERG B， FLYER N. Solving PDEs with radial basis functions［J］. Acta Numerica，2015，24：215-258.

［23］ FORNBERG B， MARTELJ M. On spherical harmonics based numerical quadrature over the surface of a sphere［J］. Advances in Computational Mathematics， 2014，40：1169-1184.

［24］ BALZER M， SCHLOMER T， DEUSSEN O. Capacity-constrained point distributions： a variant of Lloyds method［J］. ACM Transactions on Graphics （TOG）， 2009，28（3）：1-8.

［25］ FLYER N， FORNBERGB， BAYONAV， et al. On the role of polynomials in RBF-FD approximations-Part Ⅰ： interpolation and accuracy［J］. Journal of Computational Physics， 2016，321：21-38.

［26］ VLASIUK O， MICHAELS T， FLYER N， et al. Fast high-dimensional node generation with variable density［J］. Computers amp; Mathematics with Applications， 2018，76（7）：1739-1757.

［27］ SALTZMAN J， BRACKBILL J. Applications and generalizations of variational methods for generating adaptive meshes［J］. Applied Mathematics and Computation， 1982，10-11（1）：865-884.

（編輯 修榮榮）

收稿日期：2023-06-28

基金項(xiàng)目：國(guó)家自然科學(xué)基金項(xiàng)目（42004093）;山東省自然科學(xué)基金項(xiàng)目（ZR2020QD050）

第一作者：魏國(guó)華（1981-），男，高級(jí)工程師，博士研究生，研究方向?yàn)榫械卣鹳Y料處理。E-mail： weiguohua297.slyt@sinopec.com。

通信作者：谷丙洛（1988-），男，副教授，博士，研究方向?yàn)榈卣鸩▊鞑ゼ澳鏁r(shí)偏移成像方法。E-mail： gubingluo@foxmail.com。

文章編號(hào)：1673-5005（2024）04-0068-12"" doi：10.3969/j.issn.1673-5005.2024.04.007