基于改進的相位目標函數估算薄層厚度

王 鵬 胡向陽 魏水建

(中國石化石油勘探開發研究院，北京 100083)

1 引言

在地震波傳播過程中，地下介質的巖性或物性變化都將引起地震波屬性特征的變化，包括速度、振幅、頻率、相位等，這些變化已經成為地球物理儲層預測、流體識別的重要依據[1]。相位作為地震記錄的重要屬性之一，攜帶了豐富的地層反射信息，如地層反射結構特征等。

相位屬性在地震解釋中有較好的應用，尤其是在斷裂解釋中。張應波[2]認為地震相位表現了豐富、形象的地質現象，并能解釋波形變面積剖面無法解釋的多種沉積現象； 畢俊鳳[3]、柏冠軍等[4]通過觀察不同頻率的相位滯后分布，較好地實現了寬頻帶地震資料的斷層精細解釋和斷裂系統識別，獲得了較高的分辨率； 王慶華[5]利用分頻突出斷層響應的最優勢頻帶和相位對地層產狀變化敏感的特點，開發了分頻相位技術，克服了常規相干技術的不適應性，在識別斷層分布特征上起到了獨特的作用；姬戰懷等[6,7]在分頻剖面上利用瞬時相位識別小斷層； Kazmi等[8]提出使用連續相位譜的曲率屬性識別斷層特征的方法，克服了現有方法的局限性并提供了更高的體曲率屬性分辨率； Alam等[9]提出了一種利用相位譜對三維數據體中斷層進行自動檢測的技術； 劉道平[10]、趙子豪等[11]利用地震記錄的相位屬性進行工區地震層序標定，使砂體標定解釋更加可靠； 趙淑紅等[12]借鑒圖像處理方法，研究利用地震記錄的相位譜進行數據恢復。

在薄層識別和厚度估計中，人們一般以地震數據的振幅譜為基礎，利用異常特征(如零點周期[13,14]、調諧厚度[15]、振幅譜梯度[16]和主頻變化[17-19]等)開展研究。高靜懷等[20]、范明霏等[21]定義了廣義S變換并將其應用于薄層的地震探測，地層厚度識別能力從λ/4提高至λ/12； Barnes[22]通過分析薄層的振幅譜響應特征研究了頻率變化與薄層厚度之間的關系； Puryear等[23]將譜反演應用于檢測儲層厚度及沉積特征； Nowak等[24]提出了基于AVO響應的地層厚度定量預測方法，提高了譜分解方法對薄層的檢測能力； 田鑫等[25]利用地震屬性方法進行薄層砂體識別及有效厚度預測。

上述方法雖取得了一定的成果，但均很難規避頂底反射系數組合和地震子波振幅譜的影響。相比于振幅信息，相位對地層厚度的變化同樣較為敏感，但由于相位譜較復雜且對地層的響應機制不夠明確，相應的研究不多。Georgy等[26]將相位譜應用于薄互層尖滅特征分析，與振幅譜相比，相位譜能夠為尖滅點確定提供更準確的信息； 蔡涵鵬等[27]在地震頻帶范圍內應用瞬時相位譜構建了地層厚度估算的目標函數，將頂底反射系數組合項與地層厚度項巧妙地分離，在地層估算過程中不需考慮反射系數大小、極性或地震子波主頻的影響。該方法具有較好的理論意義，但存在兩個問題： 正切函數有奇異值，影響了目標函數的求解； 受噪聲污染的影響較為嚴重，在一定程度上限制了該方法的應用。

事實上，相位譜具有振幅譜不具備的優勢。例如，不論反射振幅強或弱，其相位角基本維持在同一量級；地震記錄振幅譜是地震子波振幅譜與反射系數振幅譜的乘積，而前者的相位譜則是后兩者相位譜之和，求和與乘積相比相對簡單。本文基于雙反射系數組合和地震子波常相位的假設，推導了新的地震記錄相位與地層厚度關系式，并用該式進行薄層厚度的估算，與已有的相位—厚度關系式相比，新關系式規避了相位正切函數存在的奇異點問題，提高了該方法的抗噪性。模型測試及實際資料處理結果表明，該方法可以有效提高利用相位估算薄層厚度的可靠性。

2 基本原理

假設地層反射是2個反射系數，頂、底反射系數分別為r1和r2，地層厚度為Δt，雙程旅行時為2Δt，反射系數所在時刻分別為t-Δt和t+Δt，地震子波頻譜為W(f)，復合地震記錄頻譜為S(f)，不考慮地層吸收、速度頻散、入射角度等因素，根據褶積模型，則有

S(f)=W(f)[r1e-i2πf(t-Δt)+r2e-i2πf(t+Δt)]

(1)

假設地震子波相位是零相位，時刻t=0，將式(1)右項改寫為實部U(f)與虛部D(f)之和，則實部和虛部分別為

U(f)=W(f)(r1+r2)cos(2πfΔt)

(2)

D(f)=W(f)(r1-r2)sin(2πfΔt)

(3)

式(2)與式(3)相比，可得相位譜φ(f)

(4)

求相位譜的正切值，得

(5)

Ω(f,τ)=tan[φ(f)]·cot(2πfτ)

(6)

τ為掃描變量，當τ=Δt時，中間函數Ω(f,τ)=R為一恒值(標準差為0)，當τ≠Δt時，中間函數Ω(f,τ)則是頻率f的函數，不再是恒值(標準差不為0)。

以中間函數的標準差構建目標函數Λ(τ)

(7)

(8)

(9)

(10)

3 仿真數據測試

3.1 奇異值影響分析

假設有一地層處于兩個半無限介質中間，雙程旅行時為5ms，頂、底反射系數極性均為正，比值為2，地震子波是50Hz主頻的Ricker子波(薄層厚度約為λ/8)，時間采樣間隔為0.5ms，頻率采樣間隔為0.5Hz。圖1a是反射系數(紅色虛線)及其合成地震記錄(藍色實線)的歸一化顯示。圖1b是合成記錄的振幅譜，由于受地震子波本身的帶寬限制，150Hz以上的振幅較弱，可用信息在0～150Hz之間，且振幅譜上基本看不到反射系數振幅譜的痕跡。圖1c是合成記錄的相位譜，在0～300Hz之間均有值顯示，尤其能夠清晰地刻畫150～300Hz之間的信息，不受弱振幅的影響，換言之，相位譜可利用信息的頻帶范圍比振幅譜寬。圖1d是合成記錄頻譜的實部(紅色虛線)和虛部(藍色實線)。

圖1 合成記錄及其頻域變換 (a)反射系數及其合成地震記錄； (b)合成記錄的振幅譜； (c)合成記錄的相位譜； (d)合成記錄頻譜

圖2a是根據式(6)估算的中間函數，不同曲線取的τ值不同?？梢钥吹剑敠又蹬c準確值不相等(紅線或紫線)時，曲線是非平穩的，其值隨頻率有一定的變化；當τ值與準確值相等(藍線)時，曲線呈一恒定值，滿足式(7)最佳極值逼近的假設。但藍線中還存在奇異值及其附近的異常點，這些值突然增大，有可能影響最佳極值的逼近。圖2b是根據式(9)估算的中間函數，不同曲線取的τ值不同(同圖2a)?？梢钥吹剑敠又蹬c準確值不相等(紅線或紫線)時，曲線值較大且波動較強；當τ值與準確值相等(藍線)時，曲線在零值附近波動，幅度明顯變小，需要說明的是，按照理論計算，曲線值應為零，但由于數值離散計算損失了一定的精度，曲線存在一定的誤差是合理的，也是數值計算過程中必須考慮的。

不同頻帶區間呈現的曲線特征不同，在計算目標函數時就存在頻帶范圍的選取，設置固定頻帶下限為0，頻帶上限從1Hz逐漸增加到300Hz，考察不同頻帶區間下目標函數的變化趨勢。圖2c是由圖2a根據式(7)提取的目標函數值(對數顯示)，發現τ值取準確值(5ms)時藍線的值并不總是最小，換言之，5ms并不是在所有頻帶區間內都是最佳逼近層厚。如果頻帶選0～50Hz，最佳層厚是6ms(紫線)；如果頻帶選0～150Hz，最佳層厚是5ms(藍線)；如果頻帶選0～225Hz，最佳層厚是6ms(紅線)。這種情況在理論上無法預見，且在離散計算過程中是無法避免的，而圖1c中的奇異點更是擴大了這種誤差。圖2d是由圖2b根據式(10)提取的目標函數值(對數顯示)，在35Hz之前，三條曲線并不能很好地區分開，但之后藍線(對應準確值5ms)與另兩條曲線明顯分開且不再交叉。對比結果表明，與式(7)相比，式(10)的最佳逼近值較準確，且較為穩定，受頻帶區間的影響較小。

圖3是頂底反射系數極性相反、比值為-4條件下的中間函數與目標函數，與圖2的情況類似，這里不作贅述。圖2和圖3表明，目標函數的逼近效果與頂底反射系數的極性組合、比值大小無關，最佳逼近值能反映薄層的厚度特征。換言之，基于相位譜的薄層厚度估計不用考慮頂底反射系數比的影響，這是很多振幅譜方法不具備的特點。但相位方法也有假設條件，即要求地震子波為零相位、零時刻正好在兩個反射系數中間、信噪比水平不能過低等。另外，在強噪情況下，目標函數容易在奇異點處取得極小值，如何實現全局最佳一致逼近也是需要進一步探索的問題。

圖2 中間函數及目標函數 (a)根據式(6)計算的中間函數； (b)根據式(9)計算的中間函數； (c)根據式(7)提取的目標函數； (d)根據式(10)提取的目標函數

圖3 中間函數及目標函數 (a)根據式(6)計算的中間函數； (b)根據式(9)計算的中間函數； (c)根據式(7)提取的目標函數；(d)根據式(10)提取的目標函數

3.2 穩定性分析

采用楔形模型進行目標函數的穩定性分析，地層厚度隨著道號的增加而增加，雙程旅行時從1ms逐漸增大到30ms，地震子波選取30Hz主頻的Ricker子波，第13道的厚度大致對應λ/4。圖4a是合成的楔形地震記錄，圖4b是提取的各道相位譜，窗口選取半徑為200ms的Gauss函數，相位信息在200Hz以上開始出現凌亂的情況。

分別利用常規目標函數(式(7))和改進目標函數(式(10))估計各道最佳逼近薄層厚度，采用的相位譜頻帶分別為0～50Hz、0～100Hz、0～150Hz、0～200Hz、0～250Hz和0～300Hz，估算結果分別如圖5a～圖5f所示。圖5a表明，在0～50Hz頻帶內，常規式無法較好地估算小于λ/8的薄層厚度；圖5b～圖5d表明，在0～100Hz、0～150Hz和0～200Hz三組區間上，常規式和改進式結果相近，且都與準確值接近；圖5e～圖5f表明，當頻帶區間超出200Hz時，凌亂失真的相位信息影響了薄層厚度的估算，出現較大的誤差，尤其是0～300Hz區間提取的各層厚度，改進式估算的結果有三個點與準確值(藍線)存在明顯的偏離，常規式估算的結果則大部分偏離準確值，這與上節討論的結果類似，常規式更依賴于頻帶區間的合理選取。

固定頻帶0～150Hz不變，對圖4a所示的合成記錄加入一定量的隨機噪聲，考察含噪情況下常規式與改進式的穩定性。圖6是在不同噪聲水平下提取的各道最佳逼近薄層厚度，圖6a～圖6d對應的噪聲水平(最大值與無噪記錄最大值的比值)分別為0.1%、0.2%、0.3%和0.5%。綜合分析發現，噪聲的加入對兩個目標函數逼近最佳值都有很大的影響，噪聲越大，影響越強；相比常規式，改進式具有相對較好的抗噪性，可能的原因是，隨機噪聲相位與合成記錄相位疊加，有增加奇異點數并影響估算結果穩定性的可能。

圖4楔形模型 (a)合成記錄; (b)相位譜

圖5 利用不同頻帶區間相位信息估算的薄層厚度結果 (a)0～50Hz； (b)0～100Hz； (c)0～150Hz； (d)0～200Hz； (e)0～250Hz； (f)0～300Hz

圖6 不同信噪比水平下估算的薄層厚度結果 (a)噪聲水平為0.1%； (b)噪聲水平為0.2%； (c)噪聲水平為0.3%； (d)噪聲水平為0.5%

4 實際資料應用

選取的實際資料來自陸上工區，圖7a是過工區內A1、A2、A3三口井的過井剖面，目的層為1560～1600ms(圖中黑線所示)，是砂泥巖薄層，沉積相對穩定。首先，對地震記錄進行插值(插值函數選三次樣條)，將時間采樣間隔從原先的2ms變為0.2ms。其次，選取Gauss函數作為窗口，窗口寬度參數取15～25ms(約1個波長)，沿同相軸提取目的層地震記錄的相位譜，如圖7b所示。用于計算目標函數(式(10))的頻帶分別選取0～50Hz、0～100Hz、0～150Hz和0～200Hz，四組最佳逼近結果取均值輸出，作為薄層厚度的初值(如圖7c藍線所示)。顯然，薄層厚度的初值存在較大的波動和奇異值(如紅圈中部分)，需要對藍線進行平滑處理。圖7c中的紅線是藍線經過3點中值濾波和10點均值濾波得到的光滑處理結果(即最終輸出的目的層厚度值)，隨著道號的增加，薄層厚度呈先增后減的趨勢。圖7c中的三個綠色“o”標記，代表井上統計的目的層信息： A1井(第141道)對應的雙程旅行時為2.68ms； A2井(第261道)對應的雙程旅行時為6.88ms； A3井(第441道)對應的雙程旅行時為1.37ms。對比紅線發現，A1井和A2井處的結果與井上信息吻合度較好(分別為2.97ms和7.32ms)，

圖7 實際疊加剖面處理結果 (a)疊加剖面; (b)目的層相位譜; (c)薄層厚度估算結果

表明基于相位譜能一定程度上獲得目的層的厚度信息；但A3井附近的處理結果誤差較大(3.18ms)，可能原因是該目的層厚度過薄，相位譜信息的提取不夠準確，導致估算結果偏離準確值。

圖8為根據相位屬性估算所得的目的層砂體厚度分布圖。分析預測結果可知，該地區砂巖薄層的厚度大致分布范圍為3～10m，目的層砂體呈連片條帶狀分布，根據已有地質認識，工區的物源方向為東南—西北向，切片上展示的砂體厚度變化及延展方向符合已知地質規律。表1為預測結果與井上信息進行誤差分析，T1、T2兩口驗證井的預測厚度基本在井上標準值上下浮動，相對誤差在可接受范圍內(除A3井外)，證明該方法在薄層預測方面應用效果具有一定的可信度和可行性。

圖8 目的層砂體厚度預測分布圖表1 砂體厚度預測誤差表

井名井上厚度/m預測厚度/m相對誤差/(%)A13.103.4410.97A26.707.075.52A32.475.73131.98T16.906.831.01T29.608.838.02

5 結論

(1)本文探討了相位與地層厚度的關系，并嘗試基于相位進行薄層厚度的估算，提出新的目標函數，改善常規目標函數存在的奇異值問題，在不同頻帶和信噪比水平的條件下，改進式能獲得比常規式相對更穩定和可信的估算結果。

(2)在基于相位屬性考察薄層厚度時，可以較好地規避頂、底反射系數組合的問題，這是振幅譜方法較難規避的難點；但相位方法也有假設條件，即要求地震子波為零相位、零時刻正好在兩個反射系數中間、信噪比水平不能過低等。另外，在隨機噪聲的影響下，目標函數容易在奇異點處取得極小值，如何實現全局最佳一致逼近、提高抗噪能力也是需要進一步探索的問題。

(3)當地層厚度較薄時，相位譜的估計受到干擾，無法較準確地提取薄層厚度，這是本文方法存在的一個問題，也是進一步研究的重點。

(4)地震記錄的相位譜是由地震子波相位譜與反射系數相位譜疊加得到的，與振幅譜信息相比，相位譜受弱反射的影響較小，其有效頻帶比振幅譜寬，所能利用的頻率信息更豐富，但地震相位的響應機制需要更多的研究和研討。

[1] 陳雙全,李向陽.應用局部頻率屬性檢測天然氣藏.石油地球物理勘探,2012,47(5):754-757. Chen Shuangquan,Li Xiangyang.Gas reservoir detection based on local frequency attributes.OGP,2012,47(5):754-757.

[2] 張應波.地震波相位信息的分析與應用.石油物探,1988,27(4):47-63. Zhang Yingbo.Analysis and application of seismic phase information.GPP,1988,27(4):47-63.

[3] 畢俊鳳.優勢頻帶分頻相位分析技術的應用.工程地球物理學報,2008,5(2):192-195. Bi Junfeng.The application of analyzing spectral decom-position phase of preponderant frequency channel.Chinese Journal of Engineering Geophysics,2008,5(2):192-195.

[4] 柏冠軍,趙汝敏,楊松嶺等.優勢頻帶相位分析技術在AA區塊斷層解釋中的應用.石油物探,2011,50(5):513-516. Bai Guanjun,Zhao Rumin,Yang Songling et al.Application of dominant frequency band phase analysis technique on the fault interpretation in AA block.GPP,2011,50(5):513-516.

[5] 王慶華.分頻相位分析技術在斷層解釋中的應用.工程地球物理學報,2013,10(6):896-899. Wang Qinghua.The application of frequency division phase analysis technology to fault interpretation.Chinese Journal of Engineering Geophysics,2013,10(6):896-899.

[6] 姬戰懷,嚴勝剛.分頻瞬時相位識別小斷層方法及應用.中國石油學會2015年物探技術研討會,2015,502-505. Ji Zhanhuai,Yan Shenggang.The application of frequency division instantaneous phase analysis to minor fault interpretation.Chinese Petroleum Society 2015 Geophysical Exploration Technology Conference,2015,502-505.

[7] 姬戰懷,嚴勝剛,張蓮葉.用分頻地震數據的瞬時相位識別小斷層.SPG/SEG北京2016國際地球物理會議,2016,162-165. Ji Zhanhuai,Yan Shenggang,Zhang Lianye.Minor fault detection based on instantaneous phase analysis of frequency division seismic data.SPG/SEG Beijing 2016 International Geophysical Conference,2016,162-165.

[8] Kazmi H,Alam A,Ahmad S et al.Fault characterization using volumetric curvature from continuous phase spectrum.SEG Technical Program Expanded Abstracts,2012,31:1-5.

[9] Alam A,Taylor J D.Dip azimuth and fault from continuous phase spectrum.SEG Technical Program Expanded Abstracts,2006,25:998-1002.

[10] 劉道平.儲層預測中準確的相位標定技術.石油物探,1997,36(增刊):27-32. Liu Daoping.Application of phase demarcation in re-servoir prediction.GPP,1997,36(S):27-32.

[11] 趙子豪,李凌,馬躍華等.辮狀河沉積儲層預測技術——以大港探區孔店油田為例.石油地球物理勘探,2017,52(1):152-159. Zhao Zihao,Li Ling,Ma Yuehua et al.Braided river sedimentary reservoir prediction:an example of Kongdian Oilfield in Dagang.OGP,2017,52(1):152-159.

[12] 趙淑紅,朱光明.相位在地震勘探中的用途及研究.石油物探,2005,27(4):87-89. Zhao Shuhong,Zhu Guangming.Application and research of phase in seismic prospecting.GPP,2005,27(4):87-89.

[13] Partyka G A,Gridley J,Lopez J.Interpretational aspects of spectral decomposition in reservoir characterization.The Leading Edge,1999,18(3):353-360.

[14] Partyka G A.Seismic thickness estimation:Three approac-hes,pros and cons.SEG Technical Program Expanded Abstracts,2001,20:503-506.

[15] Widess M B.How thin is a thin bed?Geophysics,1973,38(6):1176-1180.

[16] Sulistyoati T W A,Novitasari L,Winardhi S.Thickness estimation using gradient of spectral amplitude from spectral decomposition.SEG Technical Program Expanded Abstracts,2011,30:1923-1926.

[17] 孫魯平,鄭曉東,首皓等.薄層峰值頻率與厚度關系研究.石油地球物理勘探,2010,45(2):254-259. Sun Luping,Zheng Xiaodong,Shou Hao et al.The studies on relationship between thin layer seismic peak frequency and its thickness.OGP,2010,45(2):254-259.

[18] 王云專,郭雪豹,邢小林等.薄層峰值頻率特征分析.地球物理學進展,2013,28(5):2515-2523. Wang Yunzhuan,Guo Xuebao,Xing Xiaolin et al.Analysis of peak frequency characteristics of thin bed.Progress in Geophysics,2013,28(5):2515-2523.

[19] 柏冠軍,吳漢寧,趙夕剛等.地震資料預測薄層厚度方法研究與應用.地球物理學進展,2006,21(2):554-558. Bai Guanjun,Wu Hanning,Zhao Xigang et al.Research on prediction of thin bed thickness using seismic data and its application.Progress in Geophysics,2006,21(2):554-558.

[20] 高靜懷,陳文超,李幼銘等.廣義S變換與薄互層地震響應分析.地球物理學報,2003,46(4):526-532. Gao Jinghuai,Chen Wenchao,Li Youming et al.Ge-neralized S transform and seismic response analysis of thin interbeds.Chinese Journal of Geophysics,2003,46(4):526-532.

[21] 范明霏,吳勝和,曲晶晶等.基于廣義S變換模極大值的薄儲層刻畫新方法.石油地球物理勘探,2017,52(4):805-814. Fan Mingfei,Wu Shenghe,Qu Jingjing et al.A new interpretation method for thin beds with the maximum modulus based on the generalized S transform.OGP,2017,52(4):805-814.

[22] Barnes A E.Improving frequency domain thinbed analysis.SEG Technical Program Expanded Abstracts,2004,23:1929-1933.

[23] Puryear C,Castagna J P.Layer-thickness determination and stratigraphic interpretation using spectral inversion:Theory and application.Geophysics,2008,73(2):R37-R48.

[24] Nowak E J,Swan H W,Lane D.Quantitative thickness estimates from the spectral response of AVO measurements.Geophysics,2008,73(1):C1-C6.

[25] 田鑫,王緒本,張銘等.地震屬性方法在油田開發階段薄砂體識別中的應用——以印尼蘇門答臘盆地Gemah油田M油層為例.天然氣地球科學,2011,22(3):533-538. Tian Xin,Wang Xuben,Zhang Ming et al.Application of seismic attributes identifying thin sand body:An example of M Formation in Gemah Oilfield of Suma-trina Basin.Natural Gas Geoscience,2011,22(3):533-538.

[26] Georgy M,Viatcheslav P,Djalma M S F et al.Phase spectrum applied to pinch out zones analysis.9th International Congress of the Brazilian Geophysical Society,2005,1438-1443.

[27] 蔡涵鵬,龍浩,賀振華等.基于地震數據瞬時相位譜的地層厚度估算.天然氣地球科學,2014,25(4):574-581. Cai Hanpeng,Long Hao,He Zhenhua et al.Thickness estimates from instantaneous phase spectrum of poststack data.Natural Gas Geoscience,2014,25(4):574-581.