基于一階速度-應力方程的起伏地表最小二乘逆時偏移

摘要：最小二乘逆時偏移（least-squares reverse time migration，LSRTM）具有更高的成像分辨率、振幅保真性及均衡性等優勢。然而，目前的LSRTM算法大多基于水平地表假設，在面對復雜地形時無法很好地適應劇烈的起伏地表。基于二階常密度聲波方程的LSRTM算法忽略了密度變化對振幅的影響，很難在變密度介質中取得保真的成像結果。為此，從一階速度-應力方程出發，在曲線坐標系下推導相應的擾動方程和伴隨方程，并通過伴隨狀態法給出梯度更新公式，最終實現基于貼體網格的起伏地表LSRTM算法。模型試算驗證了算法的有效性和對復雜地表的適應性。結果表明，提出的算法能夠消除起伏地表的影響、壓制低頻噪聲、恢復高頻成分、均衡成像振幅，實現地下變密度介質的高分辨率和高保真度成像。

關鍵詞：貼體網格； 起伏地表； 一階速度-應力方程； 最小二乘逆時偏移

中圖分類號：P 631.4"" 文獻標志碼：A" 文章編號：1673-5005（2025）02-0108-10

Least-squares reverse time migration based on first-order velocity-stress wave equation with surface topography

LI Zhenchun^1，2， WEI Erxiang^1，2， HUANG Jianping^1，2， WANG Jiushuan³， YANG Jing³

（1.State Key Laboratory of Deep Oil and Gas， China University of Petroleum（East China）， Qingdao 266580， China;

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

3.Data Processing Center， GRI， BGP Inc， CNPC， Zhuozhou 072751， China）

Abstract： Least-squares reverse time migration（LSRTM） offers advantages such as higher imaging resolution， amplitude fidelity， and amplitude balancing. However， most existing LSRTM algorithms assume a horizontal surface and struggle to handle complex topographic variations. Additionally，" LSRTM based on the second-order constant-density acoustic wave equation neglects the impact of density variations on amplitude， making it challenging to achieve high-fidelity imaging in variable-density media. To address these limitations， we derive the perturbation and adjoint equations from the first-order velocity-stress equation under a curvilinear coordinate system. Using the adjoint-state method， we formulate the gradient update equation and implement an LSRTM algorithm based on a body-fitted grid to accommodate undulating surfaces. The results confirm the effectiveness of the proposed algorithm and its adaptability to complex topographies. Specifically， the algorithm successfully mitigates the effects of irregular terrain， suppresses low-frequency noise， recovers high-frequency components， and balances imaging amplitudes， thereby achieving high-resolution， high-fidelity imaging of subsurface variable-density media.

Keywords： body-fitted grid; surface topography; first-order velocity-stress equation; least-squares reverse time migration

隨著油氣勘探開發向復雜構造、復雜地表條件地區的不斷推進，劇烈的起伏地表已成為阻礙地震勘探發展的重要因素之一 ^［1-2］。常規的偏移成像算法大多是基于水平地表假設，對起伏地表的適應性不強，容易導致復雜地表條件下地下構造成像質量較差^［3-4］。因此，國內外學者對起伏地表條件下的地震波傳播過程及成像方法進行了大量研究，包括使用可變網格^［5-6］、垂向網格坐標變換^［7-9］及貼體網格^{［4，10-14］}等方法來擬合起伏地表，消除階梯狀離散形成的虛假散射波。然而，可變網格仍然屬于矩形網格，即使在起伏地表附近采用很細的網格，仍會出現階梯離散，無法從根本上消除虛假散射；垂向坐標變化在地表起伏比較劇烈的情況存在不穩定現象；而貼體網格能夠很好地擬合任意起伏地表，被廣泛地用于處理起伏地表問題。相較于其他偏移方法，逆時偏移成像使用雙程波方程進行波場延拓，具有不受地層傾角限制、對高陡構造精確成像的優勢^［15-17］。但是其偏移算子只是線性正演算子的共軛轉置，而不是它的逆。在地震數據采集不足、波場帶寬有限、地下構造復雜時，常規的逆時偏移方法只能對地下構造模糊成像，無法滿足實際工業化勘探開發中巖性勘探的需求^［18-19］。最小二乘偏移是把成像看作最小二乘框架下的反演問題，通過匹配反射波數據反演地下的反射系數模型，可有效改善地震波弱照明和深部偏移效果，得到具有高保真度和高分辨率的成像結果^［20-21］。最小二乘偏移思想早期主要被應用于射線理論^［22-23］和單程波理論^［24-25］，之后隨著計算能力的大幅提高，基于雙程波動方程的最小二乘逆時偏移（LSRTM）算法得到了逐步發展^［26-31］。然而，大多數LSRTM算法基于水平地表假設，沒有考慮起伏地表的影響，在一定程度上限制了其推廣應用^［4］；另外其理論框架大多是基于二階常密度標量聲波方程建立的，而實際地下介質是變密度的，忽略密度的變化將無法得到保真的成像結果^［32-34］。近年來，針對起伏地表的LSRTM算法得到了一定發展，主要是基于二階標量聲波方程^［4］和曲網格^［35］來實現。而一階速度-應力方程可以方便地處理非均勻變密度介質；貼體網格在起伏地表處有著很好的正交性，可以很方便地處理起伏地表問題；Lebedev網格^［3］在求解曲線坐標系下的離散波動方程時，可以避免交錯網格的插值誤差。因此筆者從一階速度-應力聲波方程出發，基于貼體網格和Lebedev網格，在曲網格坐標系下，推導相應的線性化波動方程和伴隨方程，進而利用伴隨狀態法給出梯度公式和更新流程，實現基于一階速度-應力方程的起伏地表LSRTM理論框架。最后通過模型試算驗證算法的正確性和適用性。

1 起伏地表LSRTM基本原理

1.1 曲坐標系下一階速度-應力聲波方程

在曲線坐標系下，一階速度-應力聲波方程可表示為

ut=1ρξxpξ+ηxpη，

wt=1ρξzpξ+ηzpη，

s2pt=ρξxuξ+ηxuη+ξzwξ+ηzwη+f

s.（1）

式中，u、w為質點速度；p為聲壓；ρ為介質密度；s為慢度，即速度的倒數；fs為震源；ξ、η為曲坐標下坐標。

將公式（1）寫成矩陣向量的形式：

LU=F.（2）

其中

L=

t0-1ρξxξ+ηxη

0t

-1ρξzξ+ηzη

-ρξxξ+ηxη

-ρzξzξ+ηzηs2t，

U=（u w p）T，

F=（0 0 fs）T.

式中，L為線性正演算子；U為波場；F為震源項。

由方程（1）可以看出，方程中的每個變量都在同一個網格的水平和垂直方向上存在空間偏導數，此時標準交錯網格不再滿足變量交錯分布的特點，求解不在交錯網格上的變量偏導數需要復雜的插值運算，增加了計算量，降低了模擬精度。為此，選用Lebedev網格來求解方程（1）。

1.2 起伏地表的線性化波動方程

根據地震波場的擾動理論，慢度場的平方s2可以表示為背景慢度場的平方s20和擾動慢度場的平方Δs2的疊加，即

s2=s20+Δs2. （3）

同樣的，根據地震波場的疊加理論可得^［36］，地震波原始波場U可以表示為由背景介質產生的背景波場U0和擾動介質產生的擾動波場Us的疊加，即

U=U0+Us. （4）

由于原始波場和背景波場都滿足波動方程（1），將式（3）、式（4）分別代入方程（1），與背景波場滿足的波動方程作差后應用Born近似條件，即可推導出一階聲波方程的逆時反偏移方程：

ust=1ρξxpsξ+ηxpsη，

wst=1ρξzpsξ+ηzpsη，

s20pst=ρξxusξ+ηxusη+ξzwsξ+ηzwsη-Δs20p0t. （5）

把式（5）利用矩陣向量表示為

LUs=Fs. （6）

其中

Us=（us ws ps）T，

Fs=0 0 mp0tT.

式中，Us為擾動波場。定義模型參數m為慢度平方的擾動，即m=-Δs20。從式（5）可以看出，擾動波場是背景波場與慢度擾動的相互作用作為二次震源產生的波場，與慢度擾動呈線性關系。

1.3 起伏地表LSRTM的伴隨方程和梯度方程

LSRTM是基于反演思想來獲取地下最優的介質模型，也就是讓正演模擬數據（即反偏移數據）與觀測數據殘差最小，這是一個最小范數問題，定義目標泛函為

F（m）=12

Lm-Uobs2=12Us-Uobs2. （7）

其中

Uobs=（uobs wobs pobs）T .

式中，m為真實反射模型;Uobs為地震波場的觀測記錄，代表地震波場在地下介質的實際傳播過程。

根據最小二乘理論，要使觀測數據與模擬數據殘差最小，需要計算目標泛函關于模型參數的梯度：

Fm=〈Us-Uobs，Usm〉.（8）

式中，Usm為偏導數波場，表示擾動波場相對于模型參數的敏感性。

式（8）中的偏導數波場可通過對式（6）兩邊同時對模型參數求導得到：

L" Us" m" = Fs m=G . （9）

其中

G" =0 0 p0tT.

直接求取式（9）中的偏導數波場是不可取的，通常要借助伴隨狀態法^［37］來計算目標泛函的梯度，因此本文中使用該方法推導對應的伴隨算子及相應的伴隨波場控制方程。

由伴隨狀態法可得

〈L*U*，U〉=〈U*，LU〉. （10）

式中，L*為L的伴隨算子；U*=（u* w* p*）表示伴隨波場。

將式（9）代入式（8），并應用式（10）中的關系可得

Fm=〈Us-Uobs，Usm〉=〈Us-Uobs，L^-1G〉=

〈（L^-1）*（Us-Uobs），G〉.（11）

這里定義伴隨波場為

U*=（L^-1）*（Us-Uobs）.（12）

進一步可以得到

L*U*=Us-Uobs. （13）

將L的表達式代入式（10）可得其伴隨算子L*：

L*=

-t0ρξxξ+ηxη

0-t

ρξzξ+ηzη

1ρξxξ+ηxη

1ρξzξ+ηzη-s2t.（14）

將式（14）代入式（13）可得伴隨波場的控制方程：

ρξxp*ξ+ηxp*η-u*t=us-uobs，

ρξzp*ξ+ηzp*η-w*t=ws-wobs，

1ρξxu*ξ+ηxu*η+1ρξzw*ξ+ηzw*η-

s2p*t=ps-pobs. （15）

比較式（1）和式（15），可以看出基于一階速度-應力方程推導出的伴隨方程與原始聲波方程存在較大差異，而基于二階位移聲波方程得到的伴隨方程與原始波動方程是相同的，這是兩種LSRTM算法的本質區別。

由式（11）和式（12）可以得到梯度更新公式：

g=Fm=p*p0t. （16）

利用共軛梯度法可以對模型進行迭代更新，其更新過程為

g^（k+1）=LT［Lm^（k）-U^obs］，

β^（k+1）=［g^（k+1）］T［Cg^（k+1）］［g^（k）］T［Cg^（k）］，

z^（k+1）=Cg^（k+1）+β^（k+1）zk，

α^（k+1）=［z^（k+1）］T［g^（k+1）］［Lz^（k+1）］T［Lz^（k+1）］，

m^（k+1）=m^（k）-α^（k+1）z^（k+1）. （17）

式中，g為最速下降法的梯度方向，可以通過對地震記錄殘差做逆時偏移成像得到；β為對應的更新步長；z和α分別為在最速下降方向和步長基礎上進行修正以后的共軛梯度方向和對應的更新步長；C為預條件算子，使用合理的正則化約束條件可以使最小二乘反演問題更加穩定，收斂速度更快，本文中使用震源波場照明作為預條件算子。

2 數值算例

2.1 起伏地表洼陷模型

為了驗證本文方法的有效性和適用性，首先采用一個起伏地表洼陷模型進行成像測試。模型如圖1所示，圖1（a）和（b）分別是物理域和計算域速度模型。模型網格點數為301×151，水平和垂直方向的網格間距均為10 m。觀測系統為全速度場接收，共301道接收，道間距為10 m；共激發60炮，炮間距為50 m；速度模型在變換到計算域后地表是水平的，因此炮點和檢波點均位于水平地表，首炮位置在（0 m，10 m）處。震源采用主頻為25 Hz的雷克子波，時間采樣間隔為1 ms，總的記錄時長為2.0 s。

圖2和圖3分別為正演模擬得到的600 ms波場快照和第30炮的炮記錄，從圖2可以看出，在起伏地表附近沒有出現虛假散射，波場干凈無頻散，說明本文正演算法在處理起伏地表方面的優勢。而受起伏地表的影響，炮記錄中的反射波同相軸發生了扭曲，在后續的成像過程中，如果忽略這個影響，將會影響成像精度。圖4（a）為慢度擾動，可以作為參結果。圖4（b）和（c）分別為RTM及其Laplace濾波結果（RTM結果是本文算法第1次迭代結果），可以看出在RTM中，淺層存在嚴重的低頻噪聲，深層同相軸的能量較弱，經過Laplace濾波后，低頻噪聲雖然得到有效壓制，但會對原有的振幅能量產生影響。圖4（d）和（e）分別為本文算法的第10次迭代和第30次迭代LSRTM結果，對比圖4（b）、（d）和（e）可以看出，隨著不斷更新計算，LSRTM逐步消除了成像剖面的低頻和高頻噪聲，使得整個剖面振幅更加均衡，成像分辨率也得到了較大提高。圖4（f）是矩形網格LSRTM第30次迭代結果，對比圖4（e）和（f）可以看出，在起伏地表附近，由于矩形網格的階梯離散導致了虛假散射，成像剖面的淺層出現了嚴重的偏移假象。圖4（g）是常密度LSRTM第30次迭代結果，與圖4（e）相比，二者成像剖面中的界面位置差異較小，需要進一步驗證密度對成像振幅的影響。

為了進一步說明該成像算法的優勢，在水平位置1.0 km處對圖4成像結果抽取單道波形，其中參考道為真實反射率模型（圖4（a））的單道波形。圖5是本文算法不同迭代次數成像結果的單道對比圖，對比圖5中的單道a、b、c可得，RTM結果的振幅與真實反射率模型有很大差別，經過Laplace濾波以后，振幅的能量與參考道仍然存在較大差異。

對比單道a、d、e可以看出，隨著迭代次數的增加，LSRTM結果越來越接近參考道振幅。圖6是常密度LSRTM和本文LSRTM結果的單道對比圖，可以看出常密度LSRTM的振幅與真實反射系數振幅存在較大差異，尤其在零振幅處存在較大的非零波動，而本文LSRTM的振幅與真實反射系數吻合度較高。以上結果說明本文LSRTM算法在處理起伏地表問題、提高成像分辨和振幅保真度方面的優勢。

2.2 起伏地表Marmousi模型

為了進一步驗證提出的基于一階速度-應力方程起伏地表LSRTM算法在成像方面的優勢，采用改進的Marmousi模型進行測試，如圖7所示，圖7（a）和（b） 分別為物理域和計算域速度模型。模型網格大小為369×199，水平和垂直方向的網格間距均為12.5 m。觀測系統采用全速度場接收，共有369道接收，道間距為12.5 m；共激發123炮，炮間距為37.5 m；速度模型在變換到計算域后地表是水平的，因此炮點和檢波點均位于水平地表，首炮位置在（0 m，12.5 m）處。震源采用主頻為25 Hz的雷克子波，時間采樣間隔為0.8 ms，總的記錄時長為3.0 s。

圖8為成像結果，圖8（a）為慢度擾動，可以作為參考結果。圖8（b）為RTM剖面，可以看到被大量的低頻噪聲覆蓋，幾乎看不到有效的構造信息。圖8（c）是圖8（b）Laplace濾波結果，雖然有效消除了低頻噪聲的污染，呈現出地下基本構造，但是與圖8（a）參考反射率模型相比，無法保持振幅的一致性。對比圖8（b）、（d）和（e）可以發現，隨著不斷的更新計算，LSRTM逐步減弱了成像剖面的低頻和高頻噪聲，基本消除了震源效應。圖8（f）是矩形網格LSRTM第30次迭代結果，對比圖8（e）和（f）可以看出，圖8（f）中成像剖面由于起伏地表的影響，在淺層出現了嚴重的偏移假象，淺層構造的有效信息被覆蓋，而深層構造的分辨率也較低。本文算法得到的成像剖面，振幅更加均衡，分辨率也得到很大的提高，成像精度大大提升。圖8（g）是常密度LSRTM第30次迭代結果，與圖8（e）相比，成像剖面在構造復雜的區域，成像精度較低，斷層以及斷層界面難以識別，說明密度變化對復雜模型的影響較大。圖9是常密度LSRTM和本文LSRTM的歸一化數據殘差收斂曲線對比，可以看出，本文算法考慮了介質密度變化，隨著迭代次數的增加，數據殘差快速減小，并逐漸收斂到5%以下。而常密度LSRTM算法忽略了密度變化，使得模擬數據與觀測數據之間存在一定差異，導致數據殘差收斂速度很慢，最終只能收斂至約40%。以上結果體現了本文算法在對起伏地表變密度介質進行偏移成像的優勢。

為了進一步說明本文算法在高分辨率成像方面的優勢，畫出圖8中本文算法不同迭代次數成像結果的總波數譜圖（圖10），即各道波數譜的疊加，其中參考曲線為真實反射率模型（圖8（a））的波數譜。對比可以看出，RTM結果中，低頻噪聲占絕對優勢，而Laplace濾波以后，低頻噪聲得到了有效壓制，但是改變了真實波數譜，與參考曲線存在很大差別。LSRTM能夠較好地保持各波數譜成分，且隨著迭代次數的增加，高頻分量逐步恢復，越來越接近真實波數譜，從而驗證了該算法能夠提高成像剖面的振幅保真度和分辨率。

3 結束語

針對起伏地表問題，采用貼體網格和Lebedev網格有限差分正演算法可以有效消除階梯離散導致的虛假散射。在此基礎上，基于一階速度-應力方程能夠處理變密度介質的優勢，推導了曲線坐標下相應的LSRTM的線性正演算子、伴隨算子和梯度更新公式，構建了一套起伏地表LSRTM算法。起伏地表洼陷模型和改進的起伏地表Marmousi模型試算結果表明，該成像算法相對于常規逆時偏移成像算法，消除了起伏地表的影響，考慮了介質的密度變化對成像剖面的影響，能夠壓制低頻噪聲、恢復高波數成分、均衡振幅，從而實現地下介質的高分辨率和高保真度成像。

參考文獻：

［1］ 孫建國.復雜地表條件下地球物理場數值模擬方法評述［J］.世界地質，2007，26（3）：345-362.

SUN Jianguo. Methods for numerical modeling of geophysical fields under complex topographical conditions： a critical review［J］. Global Geology， 2007，26（3）：345-362.

［2］ 李振春，肖建恩，曲英銘，等.時間域起伏自由地表正演模擬綜述［J］.地球物理學進展，2016，31（1）：300-309.

LI Zhenchun， XIAO Jianen， QU Yingming， et al. The review of irregular free-surface forward modeling in time domain［J］. Progress in Geophysics（in Chinese）， 2016，31（1）：300-309.

［3］ 李慶洋，黃建平，李振春，等.起伏地表貼體全交錯網格仿真型有限差分正演模擬［J］.石油地球物理勘探，2015，50（4）：633-642.

LI Qingyang， HUANG Jianping， LI Zhenchun， et al. Undulating surface body-fitted grid seismic modeling based on fully staggered-grid mimetic finite difference［J］. Oil Geophysical Prospecting， 2015，50（4）：633-642.

［4］ 李慶洋.T分布起伏地表最小二乘逆時偏移［J］.石油地球物理勘探，2019，54（5）：1075-1083.

LI Qingyang. Least-square reverse time migration on undulated surface based on T-distribution［J］. Oil Geophysical Prospecting， 2019，54（5）：1075-1083.

［5］ OHMINATO T， CHOUET B A. A free-surface boundary condition for including 3D topography in the finite-difference method［J］. Bulletin of the Seismological Society of America， 1997，87（2）：494-515.

［6］ HAYASHI K， BURNS D R，TOKSZ M N. Discontinuous-grid finite-difference seismic modeling including surface topography［J］. Bulletin of the Seismological Society of America， 2001，91（6）：1750-1764.

［7］ 董良國，郭曉玲，吳曉豐，等.起伏地表彈性波傳播有限差分法數值模擬［J］.天然氣工業，2007，27（10）：38-41.

DONG Liangguo， GUO Xiaoling， WU Xiaofeng， et al. Finite difference numerical simulation for elastic wave propagation in rugged topography［J］. Natural Gas Industry， 2007，27（10）：38-41.

［8］ TARRASS I， GIRAUD L， THORE P. New curvilinear scheme for elastic wave propagation in presence of curved topography［J］. Geophysical Prospecting， 2011，59（5）：889-906.

［9］ QU Y M， HUANG J P， LI Z C， et al. A hybrid grid method in an auxiliary coordinate system for irregular fluid-solid interface modelling［J］. Geophysical Journal International， 2017，208（3）：1540-1556.

［10］ ZHANG W， CHEN X F. Traction image method for irregular free surface boundaries in finite difference seismic wave simulation［J］. Geophysical Journal International， 2006，167（1）：337-353.

［11］ 蘭海強，劉佳，白志明.VTI介質起伏地表地震波場模擬［J］.地球物理學報，2011，54（8）：2072-2084.

LAN Haiqiang， LIU Jia， BAI Zhiming. Wave-field simulation in VTI media with irregular free surface［J］. Chinese Journal of Geophysics （in Chinese）， 2011，54（8）：2072-2084.

［12］ 丘磊，田鋼，石戰結，等.起伏地表條件下有限差分地震波數值模擬：基于廣義正交曲線坐標系［J］.浙江大學學報，2012，46（10）：1923-1931.

QIU Lei， TIAN Gang， SHI Zhanjie， et al. Finite-difference method for seismic wave numerical simulation in presence of topography—in generally orthogonal curvilinear coordinate system［J］. Journal of Zhejiang University， 2012，46（10）：1923-1931.

［13］ 黃建平，楊宇，李振春，等.基于M-PML邊界的Lebedev網格起伏地表正演模擬方法及穩定性分析［J］. 中國石油大學學報（自然科學版），2016，40（4）：47-56.

HUANG Jianping， YANG Yu， LI Zhenchun， et al. Lebedev grid finite-difference modeling for irregular free-surface and stability analysis based on M-PML boundary condition［J］. Journal of China University of Petroleum （Edition of Natural Science）， 2016，40（4）：47-56.

［14］ QU Y M， LI J L. Q-compensated reverse time migration in viscoacoustic media including surface topography［J］. Geophysics， 2019，84（4）：S201-S217.

［15］ LIU F， ZHANG G， MORTON S A， et al. An effective imaging condition for reverse-time migration using wavefield decomposition［J］. Geophysics， 2011，76（1）：S29-S39.

［16］ 王一博，鄭憶康，薛清峰，等.基于Hilbert變換的全波場分離逆時偏移成像［J］.地球物理學報，2016，59（11）：4200-4211.

WANG Yibo， ZHENG Yikang， XUE Qingfeng， et al. Reverse time migration with Hilbert transform based full wavefield decomposition［J］. Chinese Journal of Geophysics （in Chinese）， 2016，59（11）：4200-4211.

［17］ 慕鑫茹，黃建平，黎國龍，等.黏聲TTI介質中純qP波逆時偏移成像方法［J］.中國石油大學學報（自然科學版），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.

［18］ 李慶洋，黃建平，李振春.基于Students t分布的不依賴子波最小二乘逆時偏移［J］.地球物理學報，2017，60（12）：4790-4800.

LI Qingyang， HUANG Jianping， LI Zhenchun， et al. Source-independent least-square reverse time migration using student’s t distribution［J］. Chinese Journal of Geophysics （in Chinese）， 2017，60（12）：4790-4800.

［19］ 陳生昌，李代光，金成玫.高效的成像域最小二乘逆時偏移［J］.地球物理學報，2022，65（8）：3098-3107.

CHEN Shengchang， LI Daiguang， JIN Chengmei. An efficient least-square reverse time migration in image domain［J］. Chinese Journal of Geophysics （in Chinese）， 2022，65（8）：3098-3107.

［20］ TARANTOLA A. Linearized Inversion of seismic reflection data［J］. Geophysical Prospecting， 1984，32（6）：998-1015.

［21］ NEMETH T， WU C J， SCHUSTER G T. Least-squares migration of incomplete reflection data［J］. Geophysics， 1999，64（1）：208-221.

［22］ DAI W， WANG X， SCHUSTER G T. Least-squares migration of multisource data with a deblurring filter［J］. Geophysics， 2011，76（5）：R135-R146.

［23］ 黃建平，李振春，孔雪，等.碳酸鹽巖裂縫型儲層最小二乘偏移成像方法研究［J］.地球物理學報，2013，56（5）：1716-1725.

HUANG Jianping， LI Zhenchun， KONG Xue， et al. A study of least-square migration imaging method for fractured-type carbonate reservoir［J］. Chinese Journal of Geophysics （in Chinese）， 2013，56（5）：1716-1725.

［24］ 周華敏，陳生昌，任浩然，等.基于照明補償的單程波最小二乘偏移［J］.地球物理學報，2014，57（8）：2644-2655.

ZHOU Huamin， CHEN Shengchang， REN Haoran， et al. One-way wave equation least-square migration based on illumination compensation［J］. Chinese Journal of Geophysics （in Chinese）， 2014，57（8）：2644-2655.

［25］ ZHU F， HUANG J P， YU H. Least-squares Fourier finite-difference pre-stack depth migration for VTI media［J］. Journal of Geophysics and Engineering， 2018，15（2）：421-437.

［26］ 黃建平，曹曉莉，李振春，等.最小二乘逆時偏移在近地表高精度成像中的應用［J］.石油地球物理勘探，2014，49（1）：107-112.

HUANG Jianping， CAO Xiaoli， LI Zhenchun， et al. Least square reverse time migration in high resolution imaging of near surface［J］. Oil Geophysical Prospecting， 2014，49（1）：107-112.

［27］ DAI W， SCHUSTER G T. Plane-wave least-squares reverse-time migration［J］. Geophysics， 2013，78（4）：S165-S177.

［28］ DUTTA G， GIBOLI M， AGUT C， et al. Least-squares reverse time migration with local Radon-based preconditioning［J］. Geophysics， 2017，82（2）：S75-S84.

［29］ MU X R， HUANG J P， YANG J D， et al. Least-squares reverse time migration in TTI media using a pure qP-wave equation［J］. Geophysics， 2020，85（4）：S199-S216.

［30］ KIM S， KANG S-G， CHOI Y， et al. Efficient extended least-squares reverse time migration based on an excitation amplitude imaging condition［J］. Geophysics， 2023，89（1）：S31-S45.

［31］ 李振春，張閃閃，秦寧，等.TTI介質中的純qP波衰減補償最小二乘逆時偏移成像方法［J］.中國石油大學學報（自然科學版），2024，48（5）：24-35.

LI Zhenchun， ZHANG Shanshan， QIN Ning， et al. Pure qP-wave attention-compensated LSRTM in TTI media［J］. Journal of China University of Petroleum （Edition of Natural Science）， 2024，48（5）：24-35.

［32］ 李慶洋，黃建平，李振春，等.基于一階速度-應力方程的多震源最小二乘逆時偏移［J］.地球物理學報，2016，59（12）：4666-4676.

LI Qingyang， HUANG Jianping， LI Zhenchun， et al. Multi-source least-squares reverse time migration based on first-order velocity-stress wave equation［J］. Chinese Journal of Geophysics （in Chinese）， 2016，59（12）：4666-4676.

［33］ 李慶洋.復雜介質地震波正演模擬與最小二乘偏移研究［D］.青島：中國石油大學（華東），2017.

LI Qingyang. Study on seismic forward modeling and least-squares migration in complex media［D］. Qingdao： China University of Petroleum （East China）， 2017.

［34］ 郭旭.基于一階速度-應力方程的VTI介質最小二乘逆時偏移［D］.青島：中國石油大學（華東），2019.

GUO Xu. Least-squares reverse time migration based on first-order velocity-stress wave equation in VTI media［D］. Qingdao： China University of Petroleum （East China）， 2019.

［35］ QU Y M， GUAN Z， LI Z C. Topographic elastic least-squares reverse time migration based on vector P- and S- wave equations in the curvilinear coordinates［J］. Geophysical Prospecting， 2019，67（5）：1271-1295.

［36］ DUTTA G， SCHUSTER G T. Attenuation compensation for least-squares reverse time migration using the viscoacoustic-wave equation［J］. Geophysics， 2014，79（6）：S251-S262.

［37］ PLESSIX R E. A review of the adjoint-state method for computing the gradient of a functional with geophysical applications［J］. Geophysical Journal International， 2006，167（2）：495-503.

（編輯 修榮榮）

基金項目：“十四五”重大項目（2021QNLM020001）; 國家創新群體項目（41821002）; 國家自然科學基金項目（42074133）; 長慶油田公司復雜油田地震波成像技術的基礎理論研究（2023-10502）; 中石油BGP項目（2023DJ8307）

第一作者：李振春（1963-），男，教授，博士，博士生導師，研究方向為地震波傳播與成像。E-mail： leonli@upc.edu.cn。

通信作者：衛二祥（1992-），男，博士研究生，研究方向為地震波正演及偏移成像。E-mail： 976924258@qq.com。

引用格式：李振春，衛二祥，黃建平，等.基于一階速度-應力方程的起伏地表最小二乘逆時偏移［J］.中國石油大學學報（自然科學版），2025，49（2）：108-117.

LI Zhenchun， WEI Erxiang， HUANG Jianping， et al. Least-squares reverse time migration based on first-order velocity-stress wave equation with surface topography［J］. Journal of China University of Petroleum（Edition of Natural Science），2025，49（2）：108-117.