基于振幅譜互相關的薄層厚度計算

劉道理， 劉 軍， 楊登鋒， 魏旭旺， 劉 瓊

(中海石油(中國)有限公司 深圳分公司研究院，深圳 518000)

0 引言

隨著油氣勘探的不斷深入，巖性圈閉逐漸成為勘探開發的重要領域之一。相對于常規構造圈閉，巖性圈閉的隱蔽性強、儲層薄。薄層砂體厚度的有效識別和砂體尖滅的準確刻畫是巖性圈閉勘探的重點和難點[1-3]。

為了能夠準確可靠地識別薄層，很多專家在該領域做了大量的研究,Widess[4]首次利用楔狀體模型研究薄層反射與地層厚度的關系并給出了薄層的定義，將厚度小于入射子波八分之一主波長的地層定義為薄層;Neidell等[5]定義四分之一主波長對應的地層厚度為調諧厚度,在實際應用中受子波和噪聲影響，業界將厚度小于調諧厚度的地層定義為薄層;Kallweit[6-9]進一步研究了地層厚度與振幅之間的關系，得出了當地層厚度等于四分之一主波長時，反射波振幅達到最大，隨著地層厚度變薄，反射振幅也相應減小的結論，從而為利用地震波振幅信息識別薄層奠定了理論基礎；凌云等[10]利用振幅的調諧作用探測地層厚度小于1/4波長地質目標；李偉等[11]針對薄層砂體與厚層砂體不同的地震振幅響應特征，采用分頻段地震資料優選地震屬性，對秦皇島32-6油田砂體進行了預測。另外，不少的研究同樣發現地震波的頻譜特征與薄層厚度密切相關，即當地層厚度小于四分之一波長時，地層變薄，其能量變小，頻率有所提高[12-18]。根據這一特性，可以利用峰值頻率來計算薄層的厚度。Liu等[19]利用線性擬合方法近似預測薄層厚度與峰值頻率的定量關系；Sun等[20]在泰勒展開的基礎上推導出利用峰值頻率求解薄層厚度的解析式；Partyka等[21]提出譜分解技術，利用局部地震頻譜信息反演得到薄層反射系數，該項技術的缺點是在有限地震頻帶內難以準確拾取頻陷周期；Portniaguine等[22]提出了譜反演的方法，可直接反演出反射系數序列，在薄層預測中成功應用；Puryear等[23]對傳統譜反演方法進一步改進，將反射系數分解成奇偶分量，利用偶分量在厚度趨于零時的有效干涉提高了薄層反演的穩定性和分辨率；王貽朋等[24]基于匹配追蹤算法保證了高精度時頻分辨率的同時，建立地震信號偶分量的峰值頻率-薄層厚度關系曲線模板，實現譜分解薄層預測；張軍華等[2]將壓縮感知技術引入譜反演中反演反射系數，結果證實分辨率明顯提升。基于稀疏脈沖約束的常規反演技術在薄層預測中有著廣泛應用，特別是利用自然伽馬等其他測井曲線對聲波測井曲線進行曲線重構的擬聲波反演在薄層識別上效果比較理想；在井網密集區，地質統計學反演可用于薄層預測[2,25-26]。在無井或少井區這些常規的地震反演技術對薄砂體的刻畫分辨率有限。此外，基于多屬性模式識別技術同樣廣泛應用于薄砂體檢測中，模式識別方法有聚類分析、人工神經網絡、支持向量機、主成分分析等[27-29]。其基本原理是優選提取多種敏感地震屬性，根據樣本學習，在井約束的條件下可以實現薄層厚度預測，合理可靠的地震屬性優選和足夠多的學習樣本是該類方法應用的前提。

Yang Senlin等[30]提出了利用振幅譜互相關預測地震品質因子，基于該思路，筆者提出了一種利用振幅譜互相關函數計算薄層厚度的新方法。新方法首先構建目的層的振幅譜與子波振幅譜乘以反射系數振幅譜之間的相關系數函數。然后，利用高斯迭代方法使得相關系數函數達到最大值，此時對應的厚度就是薄層的厚度。通過模型測試和實際資料的應用證實新方法是可行、有效的。相較于峰值頻率方法，新方法能夠更精確地計算薄層厚度。

1 基本理論

當地震波在地層中傳播時，地層反射系數可以看成一個濾波器，地震子波為該濾波器的輸入。對于單層厚度很薄的泥包砂儲層，地震子波在其中傳播后產生的地震響應振幅譜A(f)可以表示為式(1)[12,31]。

(1)

其中：f為頻率；W(f)為子波振幅譜；r1和r2分別是薄層頂底反射系數；Δτ0為子波在該薄層中傳播的雙程旅行時。方便起見，定義薄層反射響應濾波器如式(2)所示。

(2)

其中：Δτ為子波在任意薄層中傳播過程中所需雙程旅行時。將薄層響應濾波器與震源子波相乘，得到式(3)。

(3)

薄層振幅譜和震源子波薄層響應振幅譜之間的互相關系數函數為式(4)。

(4)

為了求取相關系數C(Δτ)的最大值，需要求解如下方程：

(5)

圖1 振幅譜互相關法計算薄層厚度流程圖Fig.1 The flow chart of thin bed thickness calculation from spectral correlation

圖2 合成地震記錄及其振幅譜Fig.2 The synthetic data and its amplitude spectrum(a)合成地震記錄；(b)其振幅譜

式(5)是一個非線性方程，其解可以通過牛頓迭代法求取，具體迭代公式為式(6)。

(6)

其中：n表示迭代步數。初始迭代變量Δτ1的選取原則是：雙程旅行時大于“0”小于λ/4。迭代收斂條件為

|Δτn+1-Δτn|<ε

(7)

其中：ε為任意小的常數。

2 模型試算

圖3 振幅譜互相關函數系數曲線Fig.3 The amplitude spectral correlation function coefficient curves

圖4 楔狀體模型Fig.4 The wedge model

為了檢驗振幅譜互相關法估計薄層厚度的可行性和有效性，首先選擇三個不同厚度水平層狀薄層單砂體模型進行試算，三套單砂體對應厚度分別為5 m、10 m、20 m。砂體對應速度為3 500 m/s，密度為2 500 kg/m3；上下泥巖速度為3 000 m/s，密度為2 250 kg/m3。選用主頻40 Hz雷克子波，與模型褶積合成地震記錄，如圖2(a)所示；其中1、2、3道分別對應的砂體厚度為5 m、10 m和20 m。圖2(b)為每一道合成記錄所對應的振幅譜。基于該振幅譜且根據式(4)構建每一道所對應的互相關系數函數，其所對應函數曲線如圖3所示；其中藍、綠、紅曲線分別對應5 m、10 m、20 m單砂體厚度。由圖3可以看出，對于某一特定厚度的薄層，其互相關系數函數都有一個極大值，其所對應的雙程旅行時Δτ就是為所要求取的薄層厚度所對應的雙程旅行時。例如，圖3中紅色互相關函數曲線，在圖3上雙程旅行時約為11.4 ms處互相關函數出現極大值，當砂體速度為3 500 m/s，該極大值對應的地層厚度約為20 m。這里對比分析振幅譜互相關法和峰值頻率法兩種方法薄層估計的精度。利用振幅譜互相關法估計的三套薄層砂巖的厚度分別為4.99 m、9.99 m、20 m；對于峰值頻率法估計的厚度分別為5.03 m、10.63 m、24.30 m。從模型試算結果可見，振幅譜互相關法估計薄層厚度是可行有效的。相較于峰值頻率法，新方法估計薄層厚度更精確，這是由于新方法不需要泰勒展開近似。

圖5 基于楔狀體模型的合成地震記錄Fig.5 The synthetic data based on the wedge model

圖6 合成地震記錄所對應的振幅譜Fig.6 The amplitude spectrum of the synthetic data

圖7 利用振幅譜互相關法和峰值頻率法估計的薄層厚度Fig.7 The estimated thin bed thickness using the spectral correlation method and the peak frequency method

為了進一步驗證新方法的可行性，設計一個楔狀體模型，如圖4所示。同樣選用主頻40 Hz雷克子波，與該模型褶積合成地震記錄,如圖5所示。圖5中每一地震道對應振幅譜如圖6所示?；谡穹V互相關法和峰值頻率法計算薄層厚度如圖7所示。從圖7中可以發現，新方法計算的砂體厚度曲線與實際厚度曲線疊合的非常好；砂體厚度越大，峰值頻率法計算的砂體厚度偏離實際值越大，計算的誤差越大。相對于峰值頻率法，新方法計算精度明顯改善。該模型測試進一步證實新方法識別薄層是可行有效的。

圖8 研究區已鉆井巖性柱狀圖頻率法估計的薄層厚度Fig.8 The lithology histograms of drilled wells in research area(a)H1井；(b) H2井；(c) H3井

圖9 過研究區某一地震剖面頻率法估計的薄層厚度Fig.9 A seismic section in research area

圖10 S30層附近地震振幅譜頻率法估計的薄層厚度Fig.10 The amplitude spectrum nearby S30

圖11 提取的地震子波及其振幅譜Fig.11 The abstracted wavelet and its amplitude spectrum(a)提取的地震子波;(b)振幅譜

3 應用實例和效果

選取珠江口盆地惠州凹陷南部某區塊進行應用分析。研究區內，在地震剖面上無法直接分辨，厚度小于四分之一波長的地層都稱之為薄層。研究區珠江組S系列沉積體系以三角洲前緣近端或遠端沉積為主，砂體沉積厚度薄，巖性圈閉發育。S系列沉積的薄層砂體厚度可靠估計對于該區域油氣勘探和開發至關重要。例如S30層是該地區的主要目的層，已鉆井結果表明目的層巖性為粉砂巖，其上下地層巖性都為泥巖；地層組合特征為厚層泥巖夾單層薄砂巖，屬于典型的泥包砂類型，如圖8所示。其中，H1井S30層砂體厚度為20 m，H3井S30層砂體厚度為7.2 m。圖9位該地區原始地震資料某一任意地震剖面，其受噪聲干擾弱、信噪比高，同相軸連續性好；在S30層附近振幅譜如圖10所示，主頻約為40 Hz，有效頻帶寬度為10 Hz ～75 Hz，通過已鉆井分析目的層砂巖速度約為3 500 m/s，目的層砂體沉積厚度在四分之一波長以內，在地震剖面上難以直接分辨。S30層這套砂體沉積在地震剖面上表現為波谷反射，見圖9中藍色虛線；在剖面上很難直接識別其尖滅點及計算其沉積厚度。這里用振幅譜互相關法預測薄層厚度，新方法可靠識別薄層厚度的難點主要有：①砂體厚度越薄，預測難度越大，從圖6中可以看出，地層厚度越薄，振幅譜的差異越小，增加了預測難度；②需構建可靠的地震子波；③需要估算出可靠的地層速度。速度精度會影響厚度估算的精度，在無井區需通過對地震道集進行速度分析得到目的層段速度；在已鉆井區，通過聲波測井資料獲取目的層段速度。新方法具體步驟如下：

1)利用商業軟件在目的層段提取統計子波，提取的子波波形和振幅譜特征，如圖11所示。

2)假設薄層砂巖頂底反射強度相等，極性相反，得到泥包砂類型沉積的薄層反射響應濾波器。

3)將提取子波的振幅譜和薄層反射響應濾波器相乘，得到不同地層厚度下提取子波的響應振幅。

4)在三維地震數據體上對S30反射地震同相軸進行層位解釋追蹤，再沿著所解釋的層位上下15 ms開時窗，對S30層地震反射數據做短時傅里葉變換，得到S30薄層反射振幅譜。

5)利用實際S30薄層反射振幅譜和提取子波響應振幅譜構建得到關于地層厚度Δτ的相關系數函數。

6)求取互相關函數最大值對應的雙程旅行時。

7)由雙程旅行時和S30薄層砂體速度，計算出薄層砂體厚度。

圖12 利用振幅譜互相關法計算的S30砂體厚度與其等值線疊合圖Fig.12 The overlayed plots for the estimated thickness of S30 using the spectral correlation method and the thickness contour

圖13 S30砂體厚度等值線圖Fig.13 The thickness contour of S30

圖14 沿S30層提取的均方根振幅平面分布圖Fig.14 The RMS amplitude calculated around S30

圖15 沿S30層提取的峰值頻率平面分布圖Fig.15 The peak frequency calculated around S30

由圖12、圖13可以看出，結果顯示沉積砂體厚度普遍小于25 m，圖中暖色調區(偏紅色)指示沉積砂體厚度較厚，冷色調區(偏藍色)指示沉積砂體厚度較薄。在圖12砂體厚度平面圖的東邊，可以看出發育一套北北東向條帶狀較厚砂體沉積。為了進一步驗證所估計砂體厚度的可靠性，同樣沿著S30層提取了均方根振幅和峰值頻率屬性，分別如圖14和圖15所示。對比圖13、圖14和圖15可以看出，它們之間具有很好的對應性。砂體沉積厚的區域對應于強振幅、低峰值頻率；砂體沉積薄的區域對應于弱振幅、高峰值頻率。在圖14和圖15中同樣存在一套北北東向條帶狀沉積特征。

結合振幅、頻率屬性圖和厚度圖綜合分析S30砂體沉積微相，該套砂體地層主要有兩大沉積體系，西側靠近物源方向為三角洲前緣近端或遠端河口壩-席狀砂沉積體系；在三角洲前緣外的東南側由一系列北東向分布的條帶狀陸架砂脊沉積體系。三角洲前緣沉積不斷向海盆中心推進過程中，受波浪作用將三角洲前緣的粗粒沉積物搬運到遠離前三角洲區沉積，在沉積物供給充足的條件下，波浪作用會繼續搬運同期沉積的粗粒沉積物，從而在三角洲前緣外側形成一系列平行于岸線的砂脊，另外受后期北東-南西向潮流的沖刷和改造，導致這些砂脊條帶狀分明[32]。從沉積微相上綜合分析認為砂體沉積厚度分布較為合理可靠。此區域三口已鉆井H1、H2和H3井，應用新方法計算得到的三口井處S30層砂體厚度分別為H1井20.3 m、H2井16.5 m、H3井7.6 m；實鉆結果三口井處S30層砂體厚度分別為H1井20 m、H2井16 m、H3井7.2 m。對比計算砂體厚度和實鉆砂體厚度，可以發現計算結果與實鉆厚度非常接近，誤差很小。通過實際地震資料應用，證實了振幅譜互相關法計算薄層砂體厚度是可行、有效的。

4 結論

1)筆者提出了基于振幅譜互相關的薄層厚度計算新方法，模型測試和實際資料應用結果表明：新方法對單一薄層進行厚度計算是可行、有效的，且計算精度較高。

2)新方法適用于大套圍巖夾單一薄地層的厚度計算，對于薄互層沉積，估算的厚度與實際值可能存在偏差。若實際應用區已有鉆井，建議將已鉆井信息應用到薄層計算中，進一步提高薄層計算精度。