基于點擴散函數的深度域褶積正演方法

中圖分類號：P631 文獻標識碼：A DOI：10.13810/j.cnki.issn.1000-7210.20240287

Abstract：Seismic wave impedance inversion is thecore of reservoir prediction，and the construction ofanacurate seismic forward model is the basis of realizing high-resolution wave impedance inversion. However， it is difficult for the traditional convolution forward models based onone-dimensionalseismic wavelet to accurately simulate and characterize the imaging profile in the depth domain when there is a drastic lateral change in the subsurface velocity，which seriously affects the accuracy and reliability of wave impedance inversion results. For this reason，this paper proposes a forward modeling method of convolution in the depth domain of non-stationary space based on point spread function. Firstly，the accurate mapping relationship between seismic imaging profile andunderground reflection coefficient is derived according to the linear Born forward modeling theory.Then，an accurate convolution forward model is constructed with the point spread function as the seis mic wavelet in the multi-dimensional depth domain.Finaly，the efficient algorithm of point spread function based on Green's function of the ray theory can greatly improve the computational eficiency of the point spread function. The correctness and efectiveness ofthe proposed method are verified by the depth-domain convolution forward modeling tests of the simple horizontal layered model and the complex Marmousi model. Keywords： depth domain，wave impedance inversion，point spread function，convolutional model

沈蔚，岳玉波，秦寧，[J].石油地球物理勘探，2025，60（4）：951-957. SHENWei，YUE Yubo，QINNing：Depth-domain convolution forward modeling methodbased on point pread function[J]. Oil Geophysical Prospecting，2025，60（4） ：951-957.

0 引言

眾所周知，地震正演模型是地震反演方法的基礎[1]。當前針對深度域地震資料的反演方法依然借鑒時間域反演的思路[2-6]。即首先利用偏移速度場將深度域成像剖面轉換到時間域，使之滿足一維Ro-binson褶積正演模型所需的線性時不變性，然后利用成熟的一維時間域反演算法進行反演，最后再將時間域波阻抗反演結果轉換回深度域進行地震解釋和儲層預測。這種基于“深一時一深\"轉換的波阻抗反演方法雖然避開了深度域反演所面臨的深變子波問題。然而，當地下速度存在劇烈的橫向變化時，深度域成像剖面中蘊含了波場照明等因素所導致的復雜運動學和動力學特征信息，此時以“水平層狀介質、自激自收采集”為基本假設的一維褶積正演理論已經無法對其進行準確的表征和描述，將會嚴重影響深度域復雜地震資料波阻抗反演結果的精度和可靠性。

為實現深度域地震資料的直接反演，張靜等8基于地震多屬性變換的方法完成深度域高分辨率的地震反演；林伯香等9通過對層狀介質深度域子波的分析，提出利用速度轉換構建以“偽\"深度子波為基礎的時不變褶積模型。此后，Hu等[10]、Singh[11]胡中平等[12]、全紫荊[13]、孫萬元等[14]基于上述思想提出了“偽\"深度域波阻抗反演方法。該方法核心是保證層間旅行時不變情況下，將各層速度以深度域地震子波的介質速度進行調整，同時對地層的厚度做相應的調整，進而得到滿足線性時不變的“偽\"深度域地震數據體。后續的波阻抗反演過程借鑒成熟的時間域反演算法，最后再將反演結果校正到原始深度域。然而，從本質上來說，“偽\"深度域反演同基于“深—時一深\"轉換的反演方法都是在變換域中消除速度的垂向變化對地震子波形態的影響，使變換后的成像剖面滿足褶積正演模型所需的線性時不變性，如果不考慮深度域和時間域之間的重采樣誤差，兩者具有一致性。因此，“偽\"深度域反演的理論基礎依然是傳統以一維子波為基礎的Robinson褶積正演模型，其并不適用于深度域復雜成像資料的反演過程。

針對上述問題，本文提出了一種以點擴散函數（PointSpreadFunction，PSF）為深度域地震子波的精確褶積正演方法。不同于張金陵等[15提出的以一維PSF為基礎的深度域地震記錄合成算法（由于忽略了速度的橫向變化，其本質上仍是一維褶積正演方法），本文基于線性Born正演理論推導了以多維PSF為深度域子波表征的精確褶積正演模型，能夠有效采樣、模擬地震觀測系統和地下照明等因素對成像子波振幅、形態的影響。在此基礎上，本文給出了基于射線理論格林函數的PSF高效算法以及褶積正演模型的計算流程，并通過模型數據的褶積進行合成測試，進一步證明本文所提出的基于PSF的深度域褶積正演方法可以實現深度域成像剖面的準確模擬，能夠為深度域地震資料的高精度反演奠定重要的理論基礎。

1方法原理

1.1基于PSF的非平穩空間褶積正演模型構建

從線性Born正演理論出發，地震記錄 d 可以表示為線性Born正演算子 L 同地下反射系數模型 ?_m 的矩陣乘積

以式（1）為基礎，常規深度偏移剖面 m^′ 可以表示為L 的共軛轉置 L^T 與地震記錄 d 的矩陣乘積

m^′=L^Td

將式（1）代人式（2）得到地下反射系數 ?_m 與深度域成像剖面 m^′ 的精確映射關系

m^′=L^TLm

式中 L^TL 為Hessian矩陣，代表了震源子波、地震觀測系統、地下照明等因素所造成的成像模糊和扭曲效應[16]，也就是說常規深度偏移成像剖面是地下反射系數經過Hessian算子模糊濾波后的結果。

由于Hessian矩陣十分龐大，致使式（3）所示的映射關系難以實用化。考慮到Hessian矩陣的有效元素主要集中在主對角線附近，因此可以使用PSF來表征局部化的Hessian矩陣。此時，式（3）可表示為

m^′=Hm

式中 H≈L^TL 為點擴散函數，其為地下隨空間位置變化的多維非平穩空間信號，也是深度域地震子波的精確表征形式。以式（4）為基礎，便可以構建基于PSF的非平穩空間褶積正演模型，能夠將深度域成像剖面表征為地下反射系數同PSF的非平穩空間褶積，進而實現深度域成像剖面的精確模擬[17]。

1.2基于射線理論格林函數的PSF快速算法

傳統PSF計算方法需要對單位反射系數剖面進行串聯的正演 + 偏移運算[18]，不僅需要接近兩次常規偏移的計算成本，其空間采樣還存在明顯的限制。為保證本文方法的計算效率，在此使用岳玉波等[19]提出的基于射線理論格林函數的PSF快速算法。

以射線理論格林函數為基礎，線性Bom正演算子 L 可以表示為

G（x₀，x_s，ω）dx₀

式中： x_s 和 x_r 分別代表震源和接收點的空間位置；F（ω） 是震源子波函數， w 為頻率； m（x₀） 是散射點 x₀ 處的反射系數[20]； V 代表地下散射點的集合； G（x，x^′，ω） 是震源位于 x^′， 觀測點位于 x 的格林函數。通過求取式（5）的共軛轉置，可以得到相應的共軛偏移算子

G^*（x，x_s，ω）d（x_r，x_s，ω）dx_rdx_s

式中： m^′（x） 代表地下 x 點處的成像結果；符號 ? 代表復共軛運算。

將式（5）代人式（6），可以得到成像結果 m^′（x） 與地下反射系數剖面 m（x₀） 間的關系

式中 H（x，x₀） 為PSF，具有如下的表達形式

G（x₀，x_r，ω）G（x₀，x_s，ω）dx_rdx_s

上述PSF快速算法的實現過程主要分為以下兩步：首先，對地震數據進行循環，計算式（12）中的振幅項，并按照射線參數累加到相應的照明項中；其次，對不同的射線參數計算式（11）中平面波 Δx ），再乘以相應的照明項 A（x，p_mⁿ） 后進行累加，最終得到 x 處的PSF。

1.3基于PSF的褶積正演模型計算流程

以上述算法為基礎，即可構建基于PSF的褶積正演計算流程（圖1，進而實現深度域成像剖面的精確合成。首先，輸入道頭數據、模型速度場以及數據中提取的子波，利用PSF快速算法計算求取地下的PSF場；其次，將PSF場同地下反射系數場進行非平穩空間褶積，即可得到基于PSF的深度域褶積合成剖面。

圖1基于PSF的褶積正演模型計算流程圖

2模型測試

本文采用簡單的水平層狀模型和復雜的Marmousi模型進行測試，利用基于PSF的褶積正演模型進行深度域成像剖面的合成模擬，并將合成結果同常規Kirchhoff深度偏移剖面以及傳統基于“深—時一深\"轉換的一維褶積剖面進行對比，分析本文方法對深度域成像的褶積正演精度以及相對于傳統一維褶積正演方法的優勢。

2.1 水平層狀模型

模型網格點數為 369×375 ，縱橫向采樣間隔分別為 20m 和 8m ，模型的層速度由淺至深依次為 以及 3500m/s ，反射界面深度位于 1.2km.2. 0km 處（圖2）。使用主頻為 18Hz 的雷克子波進行正演模擬，最終合成了100炮記錄，每炮201道接收，炮點間距與道間距分別為 40m 和20m 。圖3為水平層狀模型的Kirchhoff深度偏移結果。由圖3可見，常規偏移對反射界面進行了正確歸位，但由于數據覆蓋次數的影響，成像剖面的中間能量較強、兩側能量較弱。

圖2水平層狀模型圖

圖3水平層狀模型Kirchhoff深度偏移結果

圖4為本文PSF快速算法得到的PSF場，圖中隨空間變化的局部帶限空間信號為PSF，也是該點精確的深度域地震子波[16]。PSF場與圖3所示剖面具有類似的能量分布特征。將PSF場與層狀模型反射系數場進行非平穩空間褶積，得到深度域合成剖面（圖5，僅耗時0.15s，與圖3所示剖面耗時幾乎相等，且同相軸形態特征、能量分布非常接近。進一步對比兩者在 x=3km 處成像道波形，由圖6可見，合成剖面波形（紅色）與真實的深度偏移波形（藍色）一致性較好。因此，本文方法可以實現簡單平層構造深度域成像剖面的準確模擬。

圖4水平層狀模型PSF場

圖5水平層狀模型深度域合成剖面

圖6Kirchhoff與本文方法的深度偏移結果在水平層狀模型 x=3km 處成像道波形對比

2.2 Marmousi模型

Marmousi模型的速度場網格點數為 369× 375，縱橫向采樣間隔分別為 20m 和 。首先，正演模擬了135炮地震記錄，炮點位于地表，初始炮點位于 _x=2.02km 處，炮間隔為 40m ，每炮200道雙邊接收，道間隔為 20m ，震源為主頻 18Hz 的雷克子波，對正演記錄進行常規Kirchhoff深度偏移處理，得到的成像剖面如圖8所示。

其次，根據上述觀測系統和震源子波，利用本文的PSF快速算法求取PSF場。由圖9可見，由于模型速度存在劇烈的橫向變化，計算得到的PSF場為幾何形態和振幅強度隨空間位置劇烈變化的非平穩信號（圖10a～圖10d），其攜帶的復雜運動學和動力學信息必然無法用簡單的一維子波來表征和描述。

將求取的PSF場與反射系數剖面（圖11）進行非平穩空間褶積，得到深度域合成剖面（圖12），僅耗時0.14s 由圖12可見，深度域合成剖面與圖8所示剖面的同相軸形態特征和剖面能量分布的一致性較好。

圖11Marmousi模型反射系數剖面

圖12Marmousi模型深度域合成剖面

圖13和圖14分別對比了Kirchhoff與本文方法的深度偏移結果在 x=3.4km，x=5.8km 處的成像道波形，由圖可見，合成剖面波形（紅色）同真實的深度偏移波形（藍色）基本一致。因此，本文方法可以實現復雜構造成像剖面的準確模擬。

為了進一步證明本文方法的優勢，求取了基于“深一時一深\"轉換的傳統一維褶積合成剖面（圖15）。由圖15可見，以“水平層狀介質、自激自收采集”為基本假設的一維褶積正演方法無法處理Marmousi

模型復雜的地下照明變化導致成像振幅不均衡的問題，其動力學特征在橫向和縱向都同圖8所示剖面存在很大的差異。由圖16、圖17可見，Kirch

hoff與傳統方法的深度偏移結果的波形特征和振幅強度都存在很大的偏差。因此，傳統一維褶積正演方法無法實現深度域復雜構造成像剖面的準確模擬。

3結論

（1）本文提出了一種基于PSF的深度域褶積正演方法，該方法可以快速求取以PSF表示的深度域地震子波，進而高效地實現深度域成像剖面的準確合成。（2）模型測試結果證明，相比于傳統的一維褶積正演算法，該方法可以顯著提升深度域成像剖面的模擬合成精度。（3）以本文方法為基礎，可以進一步發展基于PSF的深度域地震反演方法，有望實現深度域地震資料的直接反演計算，并提高復雜構造儲層預測的精度。

參考文獻

［1］撒利明，楊午陽，姚逢昌，等．地震反演技術回顧與展望[J].石油地球物理勘探，2015，50（1）：184-202.SALiming，YANGWuyang，YAOFengchang，etal.Past，present，and future of geophysical inversion[J].Oil Geophysical Prospecting，2015，50（1）：184-202.

[2] JIANG W，CHENXH，ZHANGJ，et al.Impedanceinversionofpre-stackseismicdatain thedepthdomain[J].AppliedGeophysics，2019，16（6）：427-437.

[3] 張杰，陳學華，蔣偉.深度域地震子波提取方法綜述[J].石油物探，2021，60（3）：353-365.ZHANG Jie，CHENXuehua，JIANGWei. Review ondepth-domainseismicwaveletestimation[J].Geophysi-calProspectingforPetroleum，2021，60（3）：353-365.

[4] 李國發，王艷倉，熊金良，等．地震波阻抗反演實驗分析[J].石油地球物理勘探，2010，45（6）：868-872.LIGuofa，WANG Yancang，XIONG Jinliang，et al.Experimental analysis on seismic impedance inversion[J].OilGeophysical Prospecting，2010，45（6） ：868-872.

[5] 陳懷震，印興耀，張金強，等.基于方位各向異性彈性阻抗的裂縫巖石物理參數反演方法研究[J].地球物理學報，2014，57（10）：3431-3441.CHEN Huaizhen，YIN Xingyao，ZHANG Jinqiang，et al.Seismic inversion for fracture rock physics parameters using azimuthally anisotropic elastic impedance[J].Chinese Journal ofGeophysics，2O14，57（10）：3431-3441.

[6] 彭真，許輝群，張偉．自適應時窗多道相關的疊后地震波阻抗反演方法[J].石油地球物理勘探，2024，59（4）：828-836.PENG Zhen，XUHuiqun，ZHANG Wei.Post-stackseismicimpedance inversionmethodwithmulti-channel correlation under adaptive window[J]. Oil Geo-physical Prospecting，2024，59（4） ： 828-836.

[7]ROBINSON E A. Seismic time-invariant convolu-tional model[J].Geophysics，1985，50（12）：2742-2751.

［8］張靜，楊勤林，王天琦．深度域地震反演方法探索[J].石油地球物理勘探，2010，45（增刊1）：114-116.ZHANG Jing，YANG Qinlin，WANG Tianqi. Pro-bing for seismic inversion in depth domain[J]. Oil Geo-physical Prospecting，201O，45（S1） ：114-116.

[9] 林伯香，胡中平，薛詩桂.利用變換深度域速度函數制作深度域合成地震記錄[J].石油地球物理勘探，2006，41（6）：640-643.LIN Boxiang，HU Zhongping，XUE Shigui. Usingvelocityfunctionintransformed depthdomaintomakesynthetic seismograph in depth domain[J]. Oil Geo-physical Prospecting，2006，41（6） ： 640-643.

[10] HUZ，ZHUH，LINB，etal.Wavelet analysisandconvolution in depth domain[C]. SEG Technical Program Expanded Abstracts，2007，26： 2733-2737.

[11] SINGH Y. Deterministic inversion of seismic data inthe depth domain[J]. TheLeading Edge，2Ol2，31（5）：538-545.

[12] 胡中平，林伯香，薛詩桂．深度域子波分析及褶積研究[J].石油地球物理勘探，2009，44（增刊1）：29-33.HU Zhongping，LIN Boxiang，XUE Shigui.Depthdomain wavelet analysis and convolution studies[J].Oil Geophysical Prospecting，2009，44（S1） ： 29-33.

[13] 全紫荊．深度域波阻抗反演方法的探討[D].成都：成都理工大學，2013.QUAN Zijing. The Study of P-wave Impedance Inver-sing Method in Depth Domain[D]. Chengdu： ChengduUniversityof Technology，2013.

[14] 孫萬元，劉仕友，李洋森，等．深度域高分辨率地震資料反演應用研究——以瓊東南盆地深度域反演為例[J].地球物理學進展，2017，32（4）：1643-1649.SUN Wanyuan，LIU Shiyou，LI Yangsen，et al.Re-search of the depth domain inversion of the seismicdata：a case of depth domain inversion in Qiongdong-nan basin[J].Progress in Geophysics，2O17，32（4）：1643-1649.

[15]張金陵，徐美茹，葉月明，等．利用點擴散函數的深度域地震記錄合成方法[J].石油地球物理勘探，2019，54（4） ： 875-881.ZHANG Jinling，XU Meiru，YE Yueming，et al.Seismogram synthetizing in the depth-domain withpoint spread function[J]. Oil Geophysical Prospecting，2019，54（4）：875-881.

[16] 盛燊，王華忠，吳成梁，等.深度域保真與高分辨期望成像子波[J].石油物探，2024，63（4）：755-765.SHENG Shen，WANG Huazhong，WU Chengliang，et al.Depth domain expected imaging wavelet withhigh-fidelityand high-resolution[J].Geophysical ProspectingforPetroleum，2024，63（4）：755-765.

[17] 焦峻峰，趙愛國，廉西猛，等，基于照明補償的變“子波\"模擬成像方法[J].石油地球物理勘探，2023，58（4）：883-892.JIAO Junfeng，ZHAO Aiguo，LIAN Ximeng，et al.Variable“wavelet” simulated imaging method basedon illumination compensation[J].Oil Geophysical Prospecting，2023，58（4）：883-892.

[18] MAOWJ，DUANWG，SUNCX，etal.Elasticleast-squares Gaussian beam imagingwith point spreadfunctions[J]. IEEE Geoscience and Remote SensingLetters，2022，19：1-5.

[19] 岳玉波，秦寧，楊哲，等．基于射線理論點擴散函數的成像域最小二乘偏移[J].地球物理學報，2023，66（2）：774-783.YUEYubo，QINNing，YANG Zhe，etal.Image domainleast-squaresmigrationbasedonray-based pointspread function[J]. Chinese Journal of Geophysics，2023，66（2）：774-783.

[20]HUH，LIUY，ZHENGY，etal.Least-squaresGaussian beam migration[J].Geophysics，2O16，81（3）：S87-S100.

（本文編輯：張偉）

作者簡介

沈蔚博士研究生，2002年生；2024年獲西南石油大學勘查技術與工程學士學位；現在西南石油大學地球科學與技術學院攻讀地質資源與地質工程專業博士學位，目前主要從事地震波阻抗反演方面的學習和研究。