Zoeppritz方程疊前多參數反演及密度敏感性分析

郭 強 張宏兵* 曹呈浩 韓飛龍 尚作萍

(①河海大學地球科學與工程學院，江蘇南京 210098； ②中科院地質與地球物理研究所，北京 100029；③河海大學力學與材料學院，江蘇南京 210098)

·綜合研究·

Zoeppritz方程疊前多參數反演及密度敏感性分析

郭 強①張宏兵*①曹呈浩②韓飛龍①尚作萍③

(①河海大學地球科學與工程學院，江蘇南京 210098； ②中科院地質與地球物理研究所，北京 100029；③河海大學力學與材料學院，江蘇南京 210098)

基于Zoeppritz方程或其近似公式的疊前地震反演應用廣泛，但多參數(尤其是密度參數)的同步反演結果不穩定。通過邊界保護正則化和馬爾科夫隨機域的軟約束建立反演目標函數，采用精確Zoeppritz方程進行反演，并使用快速模擬退火算法解非線性優化問題。數值結果表明，在小角度時，密度參數對角度變化不敏感，但對反射系數絕對值變化的貢獻度最大；二維合成數據測試結果表明，該反演算法采用小角道集可以反演出滿意的密度結果； 實際資料反演結果提供了詳細的地層信息，且與測井數據基本吻合。

疊前多參數反演 密度敏感性 Zeoppritz方程 邊界保護正則化 馬爾科夫隨機域

1 引言

地震反演是一種最有效的地震數據定量解釋方法之一，根據其采用的正演物理模型，主要分為褶積模型反演和全波形反演[1,2]。褶積模型反演，無論在疊前還是疊后，都已得到廣泛的商業應用；雖然基于聲波方程的全波形反演技術在地震成像上已有成功應用，但基于彈性波動方程的全波形反演尚在理論探索階段[2]。此外，各種最優化算法在地震反演中得到廣泛應用，如廣義線性反演法[3]，共軛梯度法[4]和模擬退火法[5]等。其中快速模擬退火法[6，7]憑借其不依賴初始模型和不易陷入局部最小值的優勢，已在很多領域得到成功應用。

目前基于AVO理論的疊前地震反演應用廣泛，它直接利用信息豐富的疊前道集數據，可以同步獲得縱波速度、橫波速度和密度[8]。但疊前地震反演仍然存在諸多問題，尤其是多參數反演結果的不穩定性問題[9，10]。由于精確Zoeppritz方程的表達式計算過程復雜，現有的文獻及商業化軟件多使用其近似式計算疊前道集反射系數，如Aki-Richard近似式[11]、Shuey近似式[12]和Fatti近似式[13]等，但這些近似式需要一定的假設條件，會導致計算誤差[14]， 特別不利于密度參數的準確提取[15，16]。此外，疊前地震反演需要解決反演病態問題[17]和地震帶限問題[18]，有關反演病態問題可以通過正則化方法得到解決[19-21]， 如邊界保護正則化方法[22-24]， 而帶限問題需要使用測井數據進行約束反演。

不同參數之間的差別會導致同步反演結果不穩定，體現為反演的彈性參數對于振幅和角度變化的敏感性。密度參數反演一直是業界研究的重點，Kabir等[25]的研究表明密度參數和速度參數對反射系數變化的貢獻度呈互耦關系，其中密度在小角度時貢獻度較大，而速度在大角度時貢獻度較大，因此需要廣角道集進行反演才能解耦密度和速度信息； Zong等[26，27]指出流體模量和楊氏模量中的密度項很難被準確估算，因為密度對振幅變化不敏感，因而直接計算得到的流體模量和楊氏模量并不準確； Liang等[28]提出了改進的Fatii方程，該方程突出密度優勢并提高密度參數對角度變化的敏感性，進而可以得到更準確的反演密度結果。

研究彈性參數中的密度參數非常必要，尤其是研究多參數同步反演中的密度敏感性。首先，建立邊界保護正則化的目標函數并結合快速模擬退火算法進行最優化；其次，分析Zoeppritz方程中反射系數對密度參數和入射角變化的敏感性；最后，通過合成數據和實際資料測試驗證本文提出的方法。

2 方法原理

2.1 正演模型

疊前地震反演是基于褶積模型的反演，其理論基礎是AVO理論，即描述平面波反射和透射系數相對入射角變化的Zoeppritz方程，其矩陣表達式為[11]

(1)

其中

(2)

式中：R是反射系數；T是透射系數；v是速度；ρ是密度； 下標“P”和“S”分別表示縱波和橫波， “1”和“2”分別代表入射介質和透視介質；θ1和θ2分別為縱波的入射角和透射角；φ1和φ2分別為轉換橫波的反射角和透射角。

2.2 目標函數

通常使用褶積模型表征波阻抗反演的數學模型，即

Y=W(θ)*R(θ)+N

(3)

式中：Y為地震記錄；W(θ)為角度震源子波；R(θ)為地震反射系數序列；N為隨機噪聲。一般情況下，假設N服從高斯分布，其數學期望為零，協方差為σ。考慮到疊前反演是一個病態問題，需要在目標函數中引入先驗信息約束，采用懲罰最小平方或者最大后驗概率估計，則目標函數可以表示為測量數據估算可信度項、先驗信息項和約束項之和，即

J(Z)=J1(Z)+λ1J2(Z)+λ2J3(Z)

(4)

式中：J1為數據項；J2為先驗項；J3是約束項；φ[·]為勢函數；C1和C2為鄰域內數據的點集；D(Z)為模型參數的梯度值；k是平滑階次；Z為待反演的模型參數，即縱波速度、橫波速度和密度。目標函數中有3個正則化參數：用于平衡數據項和先驗項之間的相互影響的“平滑參數”λ1和λ2，用于在檢測的不連續處調節梯度值的“刻度參數”δ。此外，疊前多參數同步反演的先驗項J2由三部分組成

J2=J2P(vP)+J2S(vS)+J2D(ρ)

(5)

式中：J2P、J2S和J2D分別為vP、vS和ρ的先驗項。為此，對于三個模型參數，其對應的δ值是不相同的。在式(4)中，目標函數的先驗項和約束項是以勢函數φ的形式表達；J3項的表達形式與J2項相同，但其僅在有測井數據的地震道處表達。

332332123321S123321233233

圖1 關于像素點S的MRF的鄰域系統一階鄰域包含標有“1”的4個點，n階鄰域包括標有比n小的所有點

2.3 最優化算法

反演問題可以看作是一個使目標函數達到極小的過程。式(4)中模型參數與測量數據之間呈現高度的非線性，線性類反演方法很難使該目標函數達到全局最優解，因此，使用快速模擬退火(FSA)算法[5，6]完成目標函數極小化。FSA不需要計算目標函數的偏導數，易于實現；并且不依賴于初始模型的選取，不易陷入局部最優解。采用依賴于溫度的似Cauchy分布產生新的擾動模型，即

Z(m+1)=Z(m)+T(m)sign(ξ-0.5)×

(6)

式中：Z(m)和Z(m+1)分別為三參數當前值和擾動后的值； [Zmin，Zmax]為三參數的取值范圍；ξ為[0，1]內的隨機數；T(m)是當前的溫度值。在每次迭代中，模型參數擾動后的值被接受的概率為

(7)

式中： ΔE是目標函數在模型參數擾動前后的能量差；h是預設定的常數。可見，每次迭代中模型參數被接受(或拒絕)的概率與目標函數能量的改變量有關。

3 密度敏感性分析

在地震反演中，密度項的獲取一直是個難題。因此，基于疊前反演的理論公式(精確Zoeppritz方程)，進行多參數同步反演中的密度參數敏感性的數值分析。

為了便于分析單個參數的敏感性，將式(1)中的多參數矩陣M改寫為

(8)

通過改變一個參數而固定其他兩種參數的方法，分析兩類模型反射系數對各參數和入射角變化的敏感性。圖2為第一類模型(正反射界面)反射系數隨入射角的變化曲線。由圖2c可以看出，在小角度范圍內，曲線接近于水平(梯度為零)，并且不同密度的曲線基本平行，而隨著入射角的增加，各曲線的梯度逐漸改變。這表明，在小角度時，密度參數對角度變化不敏感，即密度參數的改變并不會影響AVO曲線的變化趨勢，因此在常規梯度類AVO反演中，采用小角道集難以準確反演密度，而需要覆蓋近中、遠炮檢距的廣角道集才能獲取密度信息。圖2a和圖2b中曲線反映出，縱橫波速度參數也在適中和較大入射角范圍內才對入射角的變化敏感。

進一步分析圖2可以發現，各參數(vP、vS和ρ)曲線在不同入射角度范圍內的錯開程度是不同的，這表明不同參數在不同角度下對反射系數的貢獻度有所差異[25，31]。由圖2c可見，改變密度參數值后分別對應的五條曲線在小角度時錯開最明顯，意味著密度在小角度時對反射系數的貢獻度(導致縱波反射系數絕對值改變的程度)最大，而隨著入射角度的增加，貢獻度降低；而縱橫波速度對反射系數的貢獻度呈現出相反的規律。圖3是第二類模型(負反射界面)的反射系數隨入射角的變化曲線，顯示出相同的變化規律。需要指出的是，以上分析結果適用于縱橫波速度比變化緩慢的地層，這符合大多數地層的特點。

由以上分析可知，在小角度時，由于密度的改變(或擾動)會產生較大的反射系數絕對值的變化，使目標函數的能量差異變化增大，從而在模擬退火算法中，接受(或拒絕)該擾動的概率會增加，進而提高算法對密度參數的敏感性，提高密度反演結果的準確性。

圖2 第一類模型縱波反射系數對縱波速度(a)、橫波速度(b)和密度(c)的敏感性分析

圖3 第二類模型縱波反射系數對縱波速度(a)、橫波速度(b)和密度(c)的敏感性分析

4 合成數據試算

二維模型為一個斷背斜構造(圖4a～圖4c)，橫向共有61個CDP，間隔為12.5m，每道共151個采樣點，采樣間隔為2ms。每個CDP有17個角道集(角度間隔為3°，角度覆蓋范圍為0°～48°)。通過Zoeppritz方程計算反射系數，子波采用主頻分別為35、33和31Hz的帶限雷克子波(分別對應小、中和大角道集)，并基于褶積模型合成地震道集。圖4d為縱波速度初始模型(與真實模型的相關系數為0.65)，密度初始模型與縱波速度初始模型相似。圖5為CDP 31的合成角道集。將第15道和第45道作為虛擬井(式(4)中的測井約束項J3)對反演過程進行約束。

使用0°～12°，24°～36°和36°～48°三組角度范圍不同的角道集進行三參數反演，結果如圖6所示。對比這些反演結果：可見縱波速度和橫波速度的反演結果之間差異較小，24°～36°角道集反演的縱波速度結果與真實模型(圖4a)吻合較好，尤其是斷層處的特征明顯，36°～48°角度的反演結果的地層分界面處不夠清晰，這可能是由于大角道集距數據受調諧效應影響的結果；密度的反演結果比縱橫波速度的反演結果差，三組道集密度反演結果與真實模型(圖3c)的相關系數分別為0.86(0°～12°)、0.76(24°～36°)和0.73(36°～48°)，小角度道集反演的密度結果精度最高。取CDP 31(圖6中虛線所示)的三組密度反演結果與真實模型進行對比(圖7)，0°～12°角道集的密度結果與真實模型吻合最好。這表明，本文算法使用小角度道集可以取得較滿意的密度反演結果，但是，如果想同步獲得準確的三參數結果，需要覆蓋近中、遠炮檢距的角道集。

圖4 縱波速度(a)、橫波速度(b)和密度(c)的真實模型及縱波速度的初始模型(d)

圖5 CDP 31的合成角道集入射角范圍為0°～48°，角度間隔為3°

圖7 CDP 31三組角道集密度反演結果對比

5 實際數據測試

圖8為中國南部M工區的一條二維任意線疊加剖面，共有1981個CDP。每個CDP的角道集均有15道，角度范圍為3°～45°，間隔為3°，道間距為12.5m，時間采樣間隔為2ms。三口井分別位于CDP 53(井1)、CDP 601(井2)和CDP 1745(井3)，用井1和井2的測井數據進行約束，井3進行驗證。

圖8 M工區二維地震剖面

圖9 縱波速度(a)及密度(b)的初始模型

應用工區內幾個主要地層的參數建立了縱波速度和密度的初始模型(圖9)。初始模型合成的地震道集與實際地震數據的相關系數僅為0.477。

在反演中，結合了高階和低階鄰域，即當退火溫度T>0.005°時為三階，當0.005°>T>0.001°時為二階，當T<0.001°時為一階。優化方法FSA的初始溫度為0.05°，終止溫度為0.00001°，溫度衰減系數為0.9，并在反演過程中，逐漸降低λ值，而逐漸增大δ值。最優化算法和正則化參數與合成數據測試相同。

井3中的三套地層為泥巖、石灰巖(目的層)和砂泥巖，對應時間段分別為1680～1730ms、1730～1803ms和1803～1850ms；井2中也可見三套地層，其中泥巖和砂泥巖與井3中的相同，對應時間段分別為1570～1600ms和1725～1780ms，但井2中的石灰巖地層與井3中不同，這給井間約束增加了難度。圖10分別為3°～15°、21°～33°和33°～45°角道集的密度反演結果。對比圖10中的三組反演結果可見，小角道集反演的密度結果對地層細節信息(圖10a中箭頭所指)的分辨率最高。此外，由CDP 1745的三組密度反演結果與井3的密度測井曲線對比(圖11)可見，三組密度反演結果與測井曲線基本吻合，其中，入射角3°～15°的密度反演結果(圖11上)

圖10 實際地震資料不同角道集的密度反演結果

在1780～1810ms目的層段與測井曲線吻合最好，但是在目的層段上部(1680～1710ms)出現奇異值，這可能是隨機算法在反演過程中的不穩定性造成的誤差。

6 結論與討論

本文基于精確Zoeppritz方程并結合邊界保護正則化和MRF鄰域建立疊前反演目標函數，并采用快速模擬退火算法進行最優化反演，著重分析了多參數同步反演中密度參數的反演問題，合成數據測試及實際資料反演取得了較好的效果。

數值分析結果表明，密度參數在小角度時對反射系數絕對值變化的貢獻程度最大，隨著角度的增加，貢獻度降低；縱橫波速度參數對反射系數的貢獻度呈相反的規律。

圖11 CDP 1745的三組角道集密度反演結果與井3密度測井曲線的對比

合成數據測試表明，本文算法可以在小角度情況下獲得滿意的密度反演結果；在實際資料反演中，獲得了帶有構造細節的密度參數模型，并與測井資料基本吻合。但三參數的反演結果相互制約，獲得較好的密度反演結果，通常是以犧牲縱橫波速度的反演精度為代價。

[1] 楊文采.地球物理反演的理論與方法.北京：地質出版社，1997.

[2] Virieux J，Operto S.An overview of full-waveform inversion in exploration geophysics.Geophysics，2009，74(6)：WCC1-WCC26.

[3] Tarantola A，Valette B.Generalized nonlinear inverse problems solved using the least squares criterion.Review of Geophysics and Space Physics，1982，20(2)：219-232.

[4] Bae H S，Pyun S，Chung W et al.Frequency-domain acoustic-elastic coupled waveform inversion using the Gauss-Newton conjugate gradient method.Geophysi-cal Prospecting，2012，60(3)：413-432.

[5] 姚姚.地球物理非線性反演模擬退火法的改進.地球物理學報，1995，38(5)：643-650.

Yao Yao.Improvement on nonlinear geophysical inversion simulated annealing.Chinese Journal of Geophysics，1995，38(5)：643-650.

[6] 張霖斌，姚振興，紀晨等.快速模擬退火算法及應用.石油地球物理勘探，1997，32(5)：645-660.

Zhang Linbin，Yao Zhenxing，Ji Chen et al.Fast simulated annealing algorithm and its application.OGP，1997，32(5)：654-660.

[7] Zhang H，Shang Z，Yang C.A non-linear regularized constrained impedance inversion.Geophysical Prospecting，2007，55(6)：819-833.

[8] Chiappa F，Mazzotti A.Estimation of petrophysical parameters by linearized inversion of angle domain pre-stack data.Geophysical Prospecting，2009，57(3)：413-426.

[9] Down J E.Seismic Paramter Estimation from AVO Inversion [D].University of Calgary，Canada，2005.

[10] 田軍，吳國忱，宗兆云.魯棒性AVO三參數反演方法及不確定性分析.石油地球物理勘探，2013，48(3)：443-449.

Tian Jun，Wu Guochen，Zong Zhaoyun.Robust three-term AVO inversion and uncertainty analysis.OGP，2013，48(3)：443-449.

[11] Aki K，Richards P G.Quantitative Seismology：Theory and Methods.W.H.Freeman and Co.，1980.

[12] Shuey R T.A simplification of the Zoeppritz equa-tions.Geophysics，1985，50(4)：609-614.

[13] Fatti J L，Smith G C，Vail P J et al.Detection of gas in sandstone reservoirs using AVO analysis：A 3-D seismic case history using the Geostack technique.Geophysics，1994，59(9)：1362-1376.

[14] 鄧煒，印興耀，宗兆云等.基于逆算子估計的高階AVO非線性反演.石油地球物理勘探，2016，51(5)：955-964.

Deng Wei，Yin Xingyao，Zong Zhaoyun et al.High-order nonlinear AVO inversion based on estimated inverse operator.OGP，2016，51(5)：955-964.

[15] Zhi L X，Chen S Q，Li X Y.Amplitude variation with angle inversion using the exact Zoeppritz equations—Theory and methodology.Geophysics，2016，81(2)：N1-N15.

[16] 黃捍東，王彥超，郭飛等.基于佐普里茲方程的高精度疊前反演方法.石油地球物理勘探，2013，48(5)：740-749.

Huang Handong，Wang Yanchao，Guo Fei et al.High precision prestack inversion algorithm based on Zoe-ppritz equations.OGP，2013，48(5)：740-749.

[17] Chemingui N，Biondi B.Seismic data reconstruction by inversion to common offset.Geophysics，2002，67(5)：1575-1585.

[18] Ghosh S K.Limitations on impedance inversion of band-limited reflection data.Geophysics，2000，65(3)：951-957.

[19] 吉洪諾夫，阿爾先寧著；王秉忱譯.不適定問題的解法.北京：地質出版社，1979.

[20] Sen M K，Roy I G.Computation of differential seismograms and iteration adaptive regularization in pre-stack waveform inversion.Geophysics，2003，68(6)：2026-2039.

[21] 張豐麒，金之鈞，盛秀杰等.貝葉斯三參數低頻軟約束同步反演.石油地球物理勘探，2016，51(5)：965-975.

Zhang Fengqi，Jin Zhijun，Sheng Xiujie et al.Bayesian prestack three-term inversion with soft low-frequency constraint.OGP，2016，51(5)：965-975.

[22] 張宏兵，尚作萍，楊長春.波阻抗反演正則參數估計.地球物理學報，2005，48(1)：181-188.

Zhang Hongbing，Shang Zuoping，Yang Changchun.Estimation of regular parameters for the impedance inversion.Chinese Journal of Geophysics，2005，48(1)：181-188.

[23] Zhang H，Shang Z，Yang C.Adaptive reconstruction method of impedance model with absolute and relative constraints.Journal of Applied Geophysics，2009，67(2)：114-124.

[24] Yan Z，Gu H.Non-linear prestack seismic inversion with global optimization using an edge-preserving smoothing filter.Geophysical Prospecting，2013，61(4)：747-760.

[25] Kabir N，Crider R，Ramkhelawan R et al.Can hydrocarbon saturation be estimated using density contrast parameter？CSEG Recorder，2006，31(6)：31-37.

[26] Zong Z，Yin X，Wu G.AVO inversion and poroelasti-city with P- and S-wave moduli.Geophysics，2012，77(6)：N17-N24.

[27] Zong Z，Yin X，Wu G.Elastic impedance paramete-rization and inversion with Young’s modulus and Poisson’s ratio.Geophysics，2013，78(6)：N35-N42.

[28] Liang Lifeng，Zhang Hongbing，Dan Zhiwei et al.Prestack density inversion using the Fatti equation constrained by the P- and S-wave impedance and density.Applied Geophysics，2017，14(1)：133-141.

[29] Charbonnier P，Blanc-Féraud L，Aubert G.Determini-stic edge-preserving regularization in computed imaging.IEEE Transactions on Image Processing，1997，6(2)：298-311.

[30] Geman D，Yang C.Nonlinear image recovery with half-quadratic regularization.IEEE Transactions on Image Processing，1995，4(7)：932-946.

[31] Huang Handong，Zhang Ruwei，Shen Guoqiang et al.Study of prestack elastic parameter consistency inversion methods.Applied Geophysics，2011，8(4)：311-318.

*江蘇省南京市河海大學地球科學與工程學院，210098。Email：hbzhang@hhu.edu.cn

本文于2016年10月3日收到，最終修改稿于2017年8月30日收到。

本項研究受國家自然科學基金項目(41374116、41674113)和中國海洋石油總公司科研項目(CNOOC-KJ125 ZDXM 07 LTD NFGC 2014-04)聯合資助。

1000-7210(2017)06-1218-08

郭強，張宏兵，曹呈浩，韓飛龍，尚作萍.Zoeppritz方程疊前多參數反演及密度敏感性分析.石油地球物理勘探，2017,52(6):1218-1225.

P631

A

10.13810/j.cnki.issn.1000-7210.2017.06.012

(本文編輯：宜明理)

郭強 博士研究生，1989年生；2012年畢業于中國地質大學(武漢)地球物理學專業，獲學士學位；現于河海大學地球科學與工程學院攻讀博士學位，主要從事地震數據反演與最優化算法等研究。